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SECTION I: INTRODUCTION 



Historical Notes 



To many engineers practicing in industry the high speed digital computer 

is an unknown quantity. Just a few years ago an engineer's chances of having 

an introduction to computing while in college were slim indeed. 

Today virtually all the major engineering schools in this country have access 
to digital computers. Within the next five years the engineer who graduates 
without an introductory course in the use of computers will be the exception. 

The growth in the use of the digital computer has been tremendous . The 
first large scale digital computers appeared in the early 1950s as an out- 
growth of interest generated through the use of punched card calculators 
such as the IBM Card Programmed Calculator. Over 4,000 digital com- 
puters were installed in 1960, varying from small engineering machines 
to very large commercial and scientific data processing systems. 

Like all technological developments, the digital computer has its historical 
antecedents. Computing itself is one of the oldest human activities, re- 
quired by all civilizations for the conduct of business and the development 
of sciences. One of the most practical inventors of recent times, Charles 
Kettering, inventor of the auto self-starter, made the statement, "You 
cannot understand a thing unless you can give it a number. " In many ways 
civilization has always been tied to this need for quantitative or numerical 
description. 

In fact, the oldest surviving written document is a set of business records 
kept by a Sumerian of Mesopotamia 5,000 years ago. 

The development of computing methods was painfully slow. Perhaps this 
was due in part to the unfortunate choice of symbols (e.g. , Roman 
numerals) and of number bases (e.g. , the Babylonian choice of the base 
60) . Gradually the Arabic symbols and the base 10 won out as modes of 
computation, although vestiges of other systems remain important (e.g. , 
Roman numerals for date inscription, and base 12 for the inches in a 
foot measure) . 

Perhaps, too, this slow march forward in computing methods was due to 
the fact that they were not greatly needed, for certainly the first big 
breakthrough occurred when the science of astronomy faced tremendous 
computing efforts in its advance. In the 16th century these efforts were 
greatly eased by the development of logarithms by Napier and Briggs . 
Considering the size of numbers used in astronomy, one can easily see why 
the logarithmic transformation, which allows multiplication and division 
to be handled as addition or subtraction, was a significant contribution. 

In 1642 Blaise Pascal, a French mathematician, built the first operating 
digital computing machine. While this device demonstrated the feasibility 
of mechanical computation, the use of such devices remained very limited. 



A most interesting personality in the science of machine computation was 
the misunderstood genius Charles Babbage of England, whose major work 
was done in the middle 1800s. During his lifetime, Babbage conceived 
and designed on paper what he termed the analytical machine . In virtually 
all respects this machine was to perform like the modern digital computer, 
except that it was to be mechanical rather than electronic . The fasci- 
nation which Babbage' s work holds for present-day computer people can 
be explained when some of his ideas are examined. 

For example, his analytical machine was to have the ability to "remember" 
1,000 numbers up to 50 digits in length. This capacity is greater than 
that of some electronic computers now in existence. His machine was to 
be controlled by punched cards, an idea which he borrowed from the use 
of cards in controlling looms. In the light of the present day role of the 
punched card in data processing, Charles Babbage seems a prophet. 

It is a certainty that if Babbage had been content to build on a smaller 
scale rather than design on such a grandiose level, the science of com- 
puting would have profited greatly in the 100 years between his analytical 
machine and the present day. As it was, however, Babbage never com- 
pleted his machine , and in general we are aware of his genius primarily 
through the writings of an admirer, one Lady Lovelace, who wrote most 
descriptively of the design and philosophy of the analytical machine. 

In truth one of her statements might stand as a preface to the study of the 
modern electronic computer: "The Analytical Engine has no pretensions 
whatever to originate anything; it can do whatever we know how to order 
it to perform. " 

Charles Babbage did successfully build what is termed a difference machine 
for the computation of tables of values. Many of the early efforts in machine 
computation took this line of construction of machines to do particular 
computing jobs. A simple example will illustrate a principle made use 
of in such machines: 

To construct a table of values for 

F = M 3 + M 2 

where M is to take on the values 0, 1, 2, 3, , set down the first 

few entries and take successive differences in F: 



M 


F 








1 


2 


2 


12 


3 


36 


4 


80 


5 


150 


6 


252 



Differences 

1st 2nd 3rd 



2 


8 




10 


14 


6 


24 


20 


6 


44 


26 


6 


70 


32 


6 


102 







The third difference is a constant. The table can be completed by working 
back from the third difference to compute further entries of F by simple 
additions . It is relatively simple to build a machine to perform these 
simple additions to successive columns. 

The principle used here has been to take successive derivatives of the 
function F. 

dF 9 

^=3M 2 + 2M 

d 2 F 

iL ^ = 6M + 2 
dM 2 

d 3 F _ 6 
dM 3 



The third derivative is the constant 6. 

A great many of the functions for which tables are desired in the fields of 
navigation, astronomy and surveying can be handled by such techniques. 
The literature of computing abounds with the application of difference 
procedures . 

However, modern computers are a far cry from the difference machines 
of Babbage's day. 

At about the same time Babbage was laboring with his "Analytical Engine," 
George Boole, an English mathematician, was laying the foundation for 
later developments with his study of the algebra of logic. Boolean algebra 
is the foundation of work in machine logic (so necessary in the design of 
digital computers) and has contributed in the logic of business contracts 
and law (equally necessary in many areas of data processing). 

In the 1930s the work of modern mathematicians, in particular the late 
John von Neumann, spurred interest in the possibilities of machine com- 
putation. Various computers — both special-purpose analog devices and 
true general-purpose electromechanical digital computers — were built 



at universities. Harvard and MIT were the centers of most of this activity. 
World War II was a catalyst to machine computing in that the military 
technological advancements demanded larger and faster computational 
facilities. By the end of the war the logic for the construction of a 
general -purpose digital computer was fairly well specified. 

As makers of computing devices for accounting and business, manufac- 
turers of business machines became involved in this effort of machine 
computation. The first large scale general -purpose digital computer was 
built by IBM for use in scientific research. Known as the IBM Automatic 
Sequence Controlled Calculator, it was presented to Harvard in August 
1944. Its capacities were not surpassed until the IBM Selective Sequence 
Electronic Calculator was dedicated in January 1948. 

The SSEC was dismantled to make room for its successor, the IBM 701, 
announced in March 1953. At this time the general belief was that a few 
large-scale machines could handle all the computing needs of the country. 
The demand for more and more machines when the first few became 
available soon made apparent the vast needs for computation in both tech- 
nology and business. Today these needs have multiplied through recog- 
nition of extensions to existing needs and the discovery of new methods 
of solution of old problems. 

How does this affect the engineer ? 

The use of computers in the field of engineering analysis has experienced 
tremendous growth in the last decade. The advent of the electronic com- 
puter has created a new approach to the solution of engineering problems . 
The costly "build and try" method is now often replaced by "construct 
mathematical model and simulate. " 

For the engineer, in some respects this represents a return to the text- 
book. He must describe the physical problem by means of a mathematical 
model. Then the computer can operate on the model and describe results. 
The real power of the computer is that it will quickly describe the operating 
results for many sets of operating conditions. 

Many articles have been printed in professional publications on the solu- 
tion of difficult engineering problems by the use of electronic computers , 
and yet a lack of knowledge and experience prevents the full use of the 
computer as an engineering tool. The average engineer is simply not 
aware of what a computer can do for him and expects either too little or 
too much from this mechanism. It is highly desirable, therefore, that 
every engineer understand how problems can be described for handling 
by electronic computers. 

It is convenient to think of computing problems as falling into one of two 
classes: 

1. Straightforward computations 

2. Iterative problems 



Straightforward computations are those which have in the past been handled 
by slide rule or paper and pencil. They also include those for which the tech- 
nique of solution, while known, has required too much computation to tackle. 

Iterative problems comprise a large area, mostly unexplored. To explain 
this, the question should be asked, "What is the function of the engineer 
today?" 

The engineer is a person who builds physical systems to do particular 
jobs. If a particular construction does not work, modifications are made 
until it does work. These modifications are often the result of vague 
ideas . Engineering is done largely by trial and error . 

A very obvious conclusion which one reaches in talking with engineers is 
that they employ a minimum of their textbook analysis techniques . Math- 
ematical approaches soon pick up dust on the engineer's worktable. Why? 
One reason is that real life simply does not fit the pat equations of the 
textbooks . Often the basic difference is that in a textbook it always seems 
possible to apply mathematical techniques and arrive at solutions of the 
kind 

x = f (y, z) 

where the value of the variable of interest can be determined simply by 
plugging in known values for y and z on the right-hand side. However, in 
real-life engineering, often the problem becomes stated as 

x = f (x, y, z) 

where x is on both sides of the equation, x cannot be solved for explicitly; 
that is, there is not enough information to remove x from the right-hand 
side. The problem is too complex. 

This problem can usually be handled mathematically by estimating the 
value of x, testing the estimate, and, if it is wrong, making a better 
estimate. Essentially, numerical analysis (in computer mathematics) is 
the science of making progressively better estimates. To estimate 
repeatedly is to "iterate." This iterative idea is the heart of the computer 
approach, whether for the simplest cam problem requiring algebra, or 
for the most complex nuclear problem involving sets of differential equa- 
tions . 

This manual describes typical problems which have been handled by engi- 
neers making use of computers . No attempt is made here to give all the 
mathematical techniques required in computer analysis, nor to explain 
fully the programming techniques required for control of electronic com- 
puters . It is felt that an introduction to actual computer problems , with 
an attempt toward classification as to type, will give the practicing 
engineers insight into what computing can do and perhaps suggest possible 
uses they can personally make of engineering computing analysis. 



The Digital Computer 



The layman's image of the electronic computer contains a pinch of mystery, 
a dash of magic and a number of false analogies to human beings. These 
misconceptions have in some cases been nurtured by computer manufac- 
turers through their use of such terminology as 

1. Memory — for that part of the computer which stores information. 

2. Nerve Center — for that part of the computer which controls its 
operation. 

3 . Brain — for that part of the computer which accomplishes arithmetic 
and logical operations. 

While these analogies do help in understanding the functions of some of the 
components of a computer, they also tend to develop a camouflage of 
complexity . 

The computer, like the automobile, airplane, nuclear reactor or any other 
engineering achievement, is the product of a group of engineers. (It should 
not be surprising that in developing this product, engineers have used other 
computers to solve some of their design problems.) As in the case of the 
automobile and airplane, a person need not be capable of designing and 
building a computer in order to use one. Few persons who drive an auto- 
mobile are capable of describing the complex of gears, rods and shafts 
that allow a rotation of the steering wheel to change the direction of travel . 
This, however, does not prevent them from becoming excellent drivers. 
Similarly, a general knowledge of the functional components of a computer 
is all that is needed to use one effectively. Once the basic concepts are 
understood, greater familiarity is easily developed through experience. 

What are these basic concepts? To answer this, consider the basic com- 
ponents of a computer. 

STORAGE 

That part of a computer which most differentiates it from a calculator 
(desk-type, slide rule, adding machine) is its storage — or, to use the 
unapproved analogy, its "memory." 



STORAGE 



A computer can store numbers, alphabetic characters and some special 
symbols. In any calculation, it is necessary to store the numbers which 
begin the calculations, intermediate results and — at least temporarily — 
final results. A desk calculator has a very limited ability to store data. 
The operator usually must enter each number as it is needed. On the 
other hand, a small computer (such as the IBM 1620) can store 20,000 



digits and call them out of storage whenever they are needed in a calcula- 
tion. 

THE STORED PROGRAM 

The ability to store large amounts of data that is to be used in or has been 
developed by calculations is only part of the function of storage. Equally 
significant is the fact that the computer also contains within storage the 
program for the calculations to be performed. 



Data 

STORAGE 

Program 



A program is made up of instructions which, when acted upon by the com- 
puter, cause it to go through the sequence of operations necessary to arrive 
at a meaningful result. A single instruction may: 

1. Cause data to be brought into storage from some external source, 
such as a card reader. 

2. Cause a specified arithmetic operation to be performed on selected 
numbers. 

3. Constitute a logical test to determine what part of the program should 
be performed next. 

4. Cause results to be sent from storage to a recording device, such as 
a typewriter. 

After both the data to be operated upon and the program which describes 
the operations are in storage, the computer is free to proceed with a 
series of calculations at its own natural speed. For present-day machines, 
this may be from 50 to 500,000 additions per second. 

This leads to another salient fact: A computer can change or modify its 
own program ! Since the program is stored in much the same form as 
data, one portion of a program may be a sequence of operations which 
will examine another portion of the same program and change it by arith- 
metic or logical manipulations . What advantage is to be gained from this ? 

First, it means that one set of instructions can be used to operate on a 
number of sets of data stored in different locations in storage . As each 
set of data is processed, the program is changed to refer to the next set 
of data. 

Second, it allows the program to be changed on the basis of conditions 
which arise during the calculations. For example, suppose a calculation 
involves the evaluation of 



y = x 2 , for x $ 1 

y = -x 2 +4x-2, for x > 1. 

At the point in the calculations where the value of x becomes greater than 
1 , that portion of the program involved in the evaluation of y can be changed 
to use the second of the two equations . This change might be accomplished 
by replacing one set of instructions with a new set stored elsewhere for 
that purpose, or it could also be done by causing the computer to select 
one of two instruction sequences on the basis of whether x is greater than 
1 or not greater than 1. 

ADDRESSES 

To be able to select the item of data or the instruction to be operated upon 
next, the computer must be provided with a means of locating the desired 
information in storage. For this reason, storage is divided into units and 
each unit is identified by an address. Different computers are designed 
with different-size units of storage. The basic unit may vary from one 
decimal digit (as in the IBM 1620) to 36 binary digits (as in the IBM 7090). 
In other computers, units of one alphameric character (that is, a decimal 
digit, letter of the alphabet, or special character) or ten decimal digits 
are used. 

Whatever the size of the basic unit, each is assigned a numerical address 
which identifies the location. Manipulation of the data stored in a location 
is accomplished through the use of the address corresponding to the loca- 
tion. 

CONTROL UNIT 

The control unit of a computer provides the means whereby the stored 
program causes the desired operations to be performed in the sequence 
specified. The control unit reads an instruction from storage, examines 
it, sets up the circuit conditions to perform the operation, and, when the 
operation is completed, reads the next instruction from storage and re- 
peats these steps. 




Generally, the first operation to be performed is manually entered into 
the control unit by the machine operator. Thereafter the action of the 
computer is completely controlled by the program in storage as inter- 
preted by the control unit. An instruction which tells the computer to 
stop can be the last one in the program. 



ARITHMETIC UNIT 

This component contains the circuitry which performs arithmetic on num- 
bers taken from storage. It usually includes a limited amount of storage 
in which to hold the operands involved in the arithmetic . 



CONTROL 




At this point it is wise to pin down the kind of arithmetic accomplished in 
a computer. The two basic types of computers — analog and digital — 
have different ways of accomplishing arithmetic. The analog computer 
performs arithmetic by measuring, the digital computer by applying a set 
of rules to the numbers involved. 

To illustrate, suppose that two numbers to be added are 56.8 and 29.5. 
A simple analog device for obtaining the sum would be a ruler divided 
into 100 units. The ruler is laid on a piece of paper and marks are made 
on the paper next to the and 56.8 points on the rule. Then the ruler is 
picked up and the zero point of the ruler is placed next to the mark on the 
paper representing the 56.8 point. Now a third mark is made on the paper 
opposite the 29.5 point on the rule. All three marks, of course, are made 
to fall on a straight line. Finally, the distance between the first mark and 
the third is measured with the rule. Because of the possibility of measure- 
ment errors inherent in such a procedure, the result obtained might per- 
haps be 86.2. 



56.8 



J I .1 I I I i i i 




X 
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29.5 
X 
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86.2 



A schoolboy will add the numbers by writing one under the other and ap- 
plying a set of rules involving such things as the sum of two digits and how 
to handle the carry when the sum is greater than 9. Thus, 

1 1 

5 6.8 

2 9.5 
8 6.3 

This is the technique used by the digital computer — that is , the digital 
computer performs arithmetic by the application of a set of rules . There 
is no error inherent in such a technique . 

BINARY ARITHMETIC 

While on the subject of arithmetic, it is well to dispense with another 
question which is bound to arise in the reader's mind: What is binary 
arithmetic? 

A problem in computer design that the engineer may have recognized by 
now is that to store and manipulate decimal digits requires some device 
capable of assuming ten stable and unique states — one state to represent 
each of the digits to 9. While such devices are available (the notched 
wheels in a desk calculator are an example) , they lack the speed of opera- 
tion and small size necessary for a computer. Furthermore, they don't 
fit into an electronic system. 

Many electronic devices are available, however, which can assume two 
stable and unique states. A tube, for example, can be conducting or non- 
conducting; a magnetic field can have clockwise or counterclockwise 
rotation; a pulse can be transmitted or its transmission can be blocked. 

The computer designer has at his disposal, therefore, a number of devices 
and techniques for operating in a number system with the base 2 — the 
binary number system. The rules for arithmetic are also much simpler 
in the binary number system than in the decimal system (base 10) which 
we are accustomed to using. 

ADDITION TABLES MULTIPLICATION TABLES 








1 
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1 
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10 








1 











1 





1 



As a result, there are two basic types of digital computers: one which 
operates entirely in the binary number system and another which codes 
decimal digits as binary numbers. For the latter, the decimal digits 
appear as their binary equivalents: 
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DECIMAL BINARY 



DECIMAL BINARY 









5 


101 


1 


1 


6 


110 


2 


10 


7 


111 


3 


11 


8 


1000 


4 


100 


9 


1001 



As far as the user is concerned, it is not necessary to be familiar with 
the binary nature of the computer. The computer is capable of translating 
the engineer's language to its own before starting a program and performs 
the reverse operation when communicating its results to the engineer. 

INPUT AND OUTPUT 

So far, the methods for getting data and programs into the computer and 
for getting the results out have been ignored. Every computer must have 
a means for communicating with the outside world. 

Typical input devices are punched card readers, paper tape readers and 
manual keyboards. Card and paper tape punches, typewriters and printers 
are typical output devices . Magnetic tapes and magnetic disks provide a 
means of storing data 



CONTROL 



I 







STORAGE 






INPUT 


OUTPUT 


w 


W 







m 



ARITHMETIC 



externally from the computer in a form which allows rapid re-entry. 

The use of input-output devices is under control of the stored program. 
For example, an instruction to read a card will cause the card reader to 
start up, feed and read one card and transmit what has been read into 
storage . 

SUMMARY 

Thus there are five basic components of a computer: 

1. A place to store data and instructions in such a way that each item can 
be located when needed. 
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Exercises 



2 . A means of controlling the overall action of the system . 

3. Devices for performing arithmetic and other operations required 
in manipulating the data to obtain desired results . 

4 . A way of getting information into the computer . 

5. A way of getting results out of the computer. 



1. The decimal system has as its base the number 



2. The numeral system commonly used is the. 

numeral system. 

3. Logarithms aid greatly in doing the arithmetic operations of 

and when working with many significant 

figures . 

4. Difference principles were used in the construction of 



5. The first large scale computers were electromechanical in design; 

the middle 1950s saw the common use of the design, 

and the computers being designed and built now have solid-state 
(e.g., transistorized) components. 

6. The general iterative problems may be stated by a formula of the 
type 

7. Three computing devices with which most engineers are familiar 
are , , and 



8. One of the difficulties in predicting the behavior of complex systems 
is that engineers have trouble maintaining the 

9. Two economic aspects of a problem which most engineers are 
concerned with are and 



10. Two types of information stored in the memory of a computer 
are and 



11. Modification of the by another set of operations 

is important for repetitive and logical abilities of a computer. 

12. The unit interrogates the memory for 

information as to what to do next. 

13. The gives the physical location of infor- 
mation in memory. 

14. An analog device computes by , a digital 

computer by applying a 
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15. The binary system is a number system of. base 



16. Because the information in memory may be changed at any time, 
the digital computer is called a 

17. The physical property of is what makes a 

device an acceptable memory component. 



13 



SECTION II: THE FORTRAN LANGUAGE 



FORTRAN (FORmula TRANslation) is a language which allows the engineer 
to communicate a problem to a computer for solution. 

FORTRAN is not the natural language of a computer, nor is it the natural 
language of the engineer. Rather, it is a compromise between the two. 
To satisfy the computer, it uses symbols that the computer can understand 
and requires that the rules for their use be closely followed. To satisfy 
the engineer, it eliminates as many of the detailed computer control opera- 
tions as possible from the job of writing programs and uses a problem 
statement format close to that of the mathematical equation. 

The engineer describes his problem in the FORTRAN language; what he 
writes is translated into the natural machine language of the computer to 
be used in obtaining the solution. The translation is accomplished by the 
computer itself with the aid of a program called the FORTRAN Processor. 
The resulting machine -language program is then ready to be used to obtain 
the solution. This diagram illustrates the procedure: 



Problem 



... in FORTRAN 
language 



COMPUTER 




FORTRAN 
Processor 



Problem 




as a machine- 
language program 



COMPUTER 



Results 



Statements are the sentences of the FORTRAN language. They may: 

1 . Define the arithmetic steps which are to be accomplished by the 
computer. 

2. Provide information for control of the computer during the execution 
of the program . 

3. Describe input and output operations necessary to bring in data and 
write the results. 
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4. Specify certain additional facts such as the dimensions of a variable 
which appears in the program with subscripts. 

Because the printer and typewriter on a computer print only in upper-case 
letters, have a limited number of different characters, and are incapable 
of properly showing exponents and subscripts, FORTRAN statements at 
first appear somewhat confusing. 

An example of an arithmetic statement as it would appear on a FORTRAN 
coding sheet is: ^ 

ROOT = ( -B+SQRTF(B**2-4.*A*C))/(2.*A) \f 

Translated, this says: 

The quantity to be known as root is equal to — that is , 
can be determined by — evaluating 



Vb2T 



-B+VB^- 4AC 
2A 

where A, B, C are given values stored within the 
computer . 

Arithmetic statements look like simple statements of equality. The right 
side of all arithmetic statements is an expression which may involve 
parentheses, operation symbols, constants, variables, and functions, 
combined in accordance with a set of rules much like that of ordinary 
algebra. The symbols + and - are employed in the usual way for addition 
and subtraction. The symbol * is used for multiplication, and the symbol 
/ is used for division. The fifth basic operation, exponentiation, is rep- 
resented by the symbol **. A**B is used to represent A to the exponent 
B (that is, A B ). 

The FORTRAN arithmetic expression 

A**B*C + D**E/F - G \S 

will be interpreted to mean 

A B C + D^_- G 
F 

That is, if parentheses are not used to specify the order of operations, 
the order is assumed to be: 

1 . exponentiation 

2. multiplication and division 

3. addition and subtraction 

Parentheses are employed in the usual way to specify order. For example 
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N 



(A(B + C)) D 

is written in FORTRAN as 

(A*(B + C)) **D 

There are just three exceptions to the ordinary rules of mathematical 
notation. These are: 

1. In ordinary notation AB means A»B or A times B. However, AB 
never means A*B in FORTRAN. The multiplication symbol cannot 
be omitted. 

2. In ordinary usage, expressions like A/B-C and A/B/C are considered 
ambiguous. However, such expressions are allowed in FORTRAN and 
are interpreted as follows: 

A/B*C means (A/B)*C 
A*B/C means (A*B)/C 
A/B/C means (A/B)/C 

Thus, for example, A/B/C*D*E/F means ((((A/B)/C)*D)*E)/F. That is, 
the order of operations is simply taken from left to right, in the same way 
that 

A+B-C+D-E 

means 

(((A + B) - C) + D) - E 

3. The expression A B is often considered meaningful. However, the 
corresponding expression using FORTRAN notation, A**B**C, is not 
allowed in the FORTRAN language. It should be written as (A**B)**C 
if (A B ) C is meant, or as A**(B**C) if A( bC ) is meant. 

Besides the ability to indicate constants (like 3.57 and 2.), simple vari- 
ables (like A and ROOT), and operations (like - and *), it is also possible 
to use functions. In the previous example, SQRTF( ) indicates the square 
root of the expression in parentheses. 

Since the number of possible functions is very large, each computing cen- 
ter will have its own list of available functions , with information about 
their use. Functions given in this list must be referred to exactly as 
indicated. 

Some functions which might appear in a typical list are: 
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-Q 



FORTRAN Symbol 

ABSF(X) | X | (absolute value of X) 

SQRTF(X) yx~ 

SINF(X) sinX 

COSF(X) cos X 

ATANF(X) arctan X 

EXPF(X) e x 

LOGEF(X) ...... log e X 

Input output statements are those used to bring data into the computer to 
be stored for processing and to send results out. Typical examples are: 

GET UST[{\frf..) R£AB~l r A,B, C This statement would cause the next card in 

the card reader to be read and the three 
numbers on it stored in locations assigned 
to the values A, B and C. 

/ ** ' -~' * ' v. ? PRINT-27-ROOT This statement would cause the number in 

storage identified as the variable ROOT to 
be printed. 

PUNCH-4 r SUMA,- SUMB~This statement would cause the two values 

SUM A and SUM B to be punched on a single 
card. 

[ I | . Similar input/output statements are included for reading and writing on 

l \r T ., ^v^.'iT'-, magnetic -tapes- and- magnetic -drums-and-far-such-operations-as- r e w i nd ing 
"""" ' (1 ' ' " or-backspacing tapes. 






The numbers which follow the statements in the examples above specify 
the format in which the input or output should appear. Usually, a few 
standard formats will be specified for use in an installation. If a pro- 
grammer wishes to use something other than a standard format (for 
example, five numbers per card, in specified columns and in a specified 
notation), a new format can be introduced by means of a format statement. 
The FORTRAN Processor interprets the format statement and causes the 
machine-language program being generated to conform to the desired 
format in handling input or output statements which refer to it. 

Control statements allow the programmer to state the flow of his program 
through use of the logical powers of the computer. In writing a FORTRAN 
program, any statement in the program which is referred to by another 
statement must be given an identifying number (every statement can be 
given a number if desired but this is not necessary). The control state- 
ments refer to these identifying numbers for the purpose of branching 
from one part of the program to another. 

Important statements in this category are illustrated by the following 
examples: 
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GO TO \ VlQ/L\ " This statement indicates that the next 

statement to be executed (after having been 
converted to a machine-language program, 
,< of course) is the one numbered 4. 

GO TO , (4V"I8"" s "~20;-40)?\ K This statement is referred to as a computed 

GO TO since the value of K is computed in 
a previous statement. If, at the time this 
GO TO were executed, K = 3, then the third 
alternate (statement 20) would be the one 
chosen. 

/ IF (A**2 - B) lG-r20T~SO~- This statement allows one of three alternate 

-I. /'""(. A-v": ' <?. <Z-Q } ~pj. (T, j { ^q -j fx \/C) » next instructions to be chosen according to 

f? J c/_i "** ■■- n ' \Ar - - ^ whether the expression in parentheses is 

w """ A •-'•' ; r" t~~ >'--'■ U 1 ' ^'^ ;■ &-0 1 O LS-OJ i ess than, equal to, or greater than zero 

£-"]..- f X' G o i V L ~£ G ( in t ' iat or< ^ er ) • 



t > 

X 



This statement is a command to execute 
the sequence of statements which follow, 
up to and including 25 (or any specified 
statement number) , repeatedly. The first 
time, let the subscript I on the variables 
appearing in the statements equal 1 (the 
first number or variable following the = 
sign). The second and each subsequent 
time the sequence is executed, increase 
I by 2 (the third number or variable following 
the equal sign — this can be left blank if it 
is to be 1). When I becomes greater than 
M (the second number, or variable as in 
this case) do not repeat, but continue by 
executing the statement after 25. 

END This is the last statement in any 'FORTRAN 

program*. fc~A_/> <* r) C' ..*.> \ ^ 

In the DO example, reference was made to the changing of subscripts on 
variables. For example, A(I, J) would refer to a two-dimensional array 
of quantities a^. More will be said about this in the next section. 

What has been said about the FORTRAN language so far will be more 
meaningful by examining a complete FORTRAN program. The complete 
FORTRAN program to find the roots of the quadratic equation 

ax 2 + bx + c = 

might be written in the following logical steps: 

Arithmetic Statements 

The two roots are evaluated as arithmetic statements. A separate state- 
ment evaluating the discriminant is used to avoid the extra computation 
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that could be involved if this expression were written in both root state- 
ments and to permit later testing of the discriminant. 

2 D =B**2 - 4^*A*C 

4 ROOT 1 = (-B+SQRTfe (D) )/ (2.-fl*A) 

5 ROOT 2 = (-B-SQRTF, (D) )/ (2;;0*A) 

All the parentheses in statements 4 and 5 are necessary. D is enclosed 
in parentheses to define the expression to be operated on by the square 
root function. -B + SQRTF (D) is in parentheses so the additions and sub- 
tractions will be performed prior to the division. 2. 0*A is in parentheses 
so the multiplication will be performed prior to the division. 

Input/Output Statements 

Statements are added to enter values for the constants A, B and C and to 
print the results of computation. Logically the read statement must pre- 
cede the arithmetic statements so that values will be available for compu- 
tation in the arithmetic expressions . The print statement must be after 
the arithmetic statements, so that values will be available for printing. 

&ET LIJTT- 

1 RRAD-Wk/ 



^(a, B, C \ 



2 D = B**2 - 4.0*A*C 

4 ROOT 1 = (-B + SQRTF (D) ) / (2. 0*A) 

5 ROOT 2 = (-B - SQRTF (D) )/ (2.0*A) 

J* UT " l|fT / \ 

6 *T>RINT-20( A, B, C, ROOT 1, ROOT 2) 

In addition to the roots, the values of the constants A, B and C are printed 
to identify the results with the input. 

Control Statements 

Statements are now added for decision making and flow of the program: 
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(S., 



LIST 



i z. 



o)tq- l? 



•1= 

STATI 



C f « COMMENT 



STATEMENT 
NUMME 



.L^ 



^A 



±A 



A 



Ai 



jj 



FORTRAN STATEMENT 



/t,£,*AJA,? 



fi,=A**A'. 



Mfahi, 



*,*AWX- 



*MTA* x < \ 



PAIATJA 



&A ?d X 



'ATATJA 



*A M ,' ■ 



¥A , , 



'''■'I 



*-A,*JAC . 



4 t ,A 



>A >\£, 



t"l I I I I I I I I 



I I I I I 



''■'■■' 



I I I I I I I I 



*ASA*.T,k<,£>,? t ),/Ct l .A*.*,) 



■*-■■■* ■ ' ' ' i \ " i i - ■ '■ - ■ ■ ■ ■ ■ ■* 

RrfGA rX (A I l/(X .A*A) 



A,A>A,A0ATX>A0ATA 



i i i i i i 



' I ' ' ' ' t c 



^,ttB,tA\ i i 



i i i i i i 



_L_I '■''■''' I I I I I L. 



I I ' ' ' ' I ' 



' ■ ' ' ■ ' ' ' 



■■' '' 



Statement 3 is added to test for the mathematical possibility of a negative 
discriminant (complex roots). If the value of D is zero or positive, the 
flow of the program will be the normal flow to statement 4. If the value 
of D is negative, the flow will be to statement 8. Several alternatives are 
possible at this point depending upon the desired results: 

1. If complex roots are desired, they may be calculated and printed in 
their real and complex parts . 

2. This particular set of data could be ignored completely by transferring 
control to statement 1, where the next set of values would be read. 

3. Statement 8 could be a STOP statement which would cause the computer 
to stop execution of the program and permit manual intervention. This 
is not a recommended procedure. 

4. Statement 8 could be a PRINT statement to signal an error condition. 
In this program values A, B and C will be printed without the root 
values and can later be analyzed to determine what was wrong with 
this set of data. 

Statements 7 and 9 are added to complete the flow of the program by re- 
turning to statement 1 to read the next set of values . 

Statement 11 is the END statement to indicate the end of this program. 

This program will continue to be repeated as long as there are cards in 
the hopper of the card reader. When an attempt is made to read a card 
after the hopper has emptied, the computer will stop. 

Additional examples of FORTRAN programs will be found in the sections 
that follow. The use of additional statements will be explained at that 
time. 
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Exercises 



The FORTRAN language is designed to simplify engineering use of IBM 
computers. The use of the language with a particular IBM computer re- 
quires the use of a FORTRAN Processor for that machine. 

The power of IBM computers increases as the size and speed, in terms of 
machine components, increase. For the more powerful machines the 
FORTRAN language is expanded to take advantage of the additional types 
of components. For example, computers with magnetic tape input/output 
units allow use of FORTRAN statements to control the use of tape . 

Should the reader wish to become more expert on the FORTRAN language 
in general, or as it applies to a particular IBM machine, he may contact 
the local IBM representative for information on available IBM manuals 
and education courses. 



Write FORTRAN-language programs to: 

1. Put the numbers 1,2, .0396, 30,000.0 and -6.5 into the storage of 
a computer. 

*2. Compute C=A+B and print the resultant value for C, where A=3.5 
and B=2.694. 

3. Compute C=A+B, where A and B are to be read into the computer and 
C is to be printed. 

*4. Compute C=A+B if B > 
C=A-B if B < 

where A and B are read into the computer, and C is to be printed. 

*5. Generate the integers 1 to 100 and print. 



6. Generate the value of sin x for x = .0000, .0001, .0002, 

.5. Print the value of x and sin x at each interval. 

*7. Change the program for solution of a quadratic equation to include 
solution for complex roots. 

* Answers to these exercises are found in Appendix I. 
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SECTION HI: SOME NEW AND OLD CONCEPTS 

Each new area of technology seems to construct a particular vocabulary 
for the necessary descriptions of efforts and communication of ideas. 
Computing is no exception to this rule; new terms such as "stored program" 
and "bit" have been added and new meaning given to old terms such as 
"iteration," "algorithm" and "memory." 

As far as use of computers is concerned, this growth in vocabulary is 
brought about by the new man-machine relationship. Control of a com- 
puting machine requires complete accuracy. This demands conciseness 
in statement of what is to be done, and preconception of what alternatives 
the computer takes as the result of all possible logical decisions it is 
asked to make. Of course, the usual give and take is allowed in regard to 
the testing of ideas. 

This section is designed to explain some of the procedures or techniques 
used in problem definition to simplify communication with the computer. 
In some cases this amounts to a re-emphasis of standard mathematical 
techniques . 

Block Diagramming 

The old cliche' — a picture is worth a thousand words — has particular 
truth for the engineer in preparing a program for a computer. Since the 
engineer normally finds many uses for pictorial representations in his 
work, the block diagram is readily accepted as a powerful tool in comput- 
er programming. 

A block diagram is a picture of the steps which must be performed to ac- 
complish a particular job. The major function of the diagram is to clarify 
what must be done as a result of each decision. All the alternatives are 
accounted for in each case . In programming a computer it must be remem- 
bered that the computer has no way of anticipating the requirements of a 
problem and must be provided with all the information needed to reach a 
solution. 

The amount of information put into the block diagram, however, will de- 
pend upon personal preference, the programming techniques used, and 
the type of computer to be used. For example, Figure 1 illustrates two 
ways of block-diagramming the solution of a quadratic equation. The de- 
tailed diagram shows every step in the computation and is the type that 
might be used if the programmer were forced to work in machine language . 
The second diagram shows only the three main steps in the program but 
is quite sufficient for a programmer using the FORTRAN language. 




To further illustrate, consider a problem which involves the repetition of 
the same operation over and over. Let the problem be to create the sum 



to 5 



//■n\ r&^i-^ 2 






This might be block-diagrammed in either of the two ways illustrated in 
Figure 2. 
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A. 
DETAILED 





[ Read j 
I A,B,C J 




V 




Compute 
X 




v 


B. 


[ Punch X ' 


GENERAL 







B' - 4AC 



D=Vb 2 - 4AC 




Figure 1. 
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ReadXj, Yj 
i=l, 2, . . , 5 



A. B. 



ReadXj, Y. 
i=l,2,..,5 



D l = ( X l " Y l) 



D 2 = (X 2 -Y 2 ) 



etc. 



D 5 = < X 5- Y 5> 



Z = Dj + D 2 
+D 3 +D 4 -fD 5 



Set i = 1, 
Z = 



Compute 
Z=Z + (X r Y.) 2 




Set i = i+1 







Punch Z 



Figure 2. 
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Subscripting 



The two diagrams illustrate an important programming concept — that of 
looping. Diagram A corresponds to a program in which each (X. - Y.) is 
computed separately and then all are summed to obtain Z. With only 
five values of i, this is not too unlikely a method of approaching the prob- 
lem. But, consider a situation in which i = 100 or 1000 or may vary 
according to some other characteristic of the problem of which this com- 
putation is a part. In such a case not only would the diagram be large and 
time-consuming to draw, but the program would also be large and time- 
consuming to write . 

In diagram B, advantage has been taken of the computer's ability to make 
logical decisions and to modify its own program . Since the number of 
times the computation is repeated depends only on the value of i (and the 
presence of successive values Xj and Y.), this one diagram — and the 
program which would be based on it — applies for any value of i. 

The algebraically incorrect statement i = i + 1 means replace the current 
value of i with i + 1. Similarly, Z = Z + (X^ - Yj)^ means replace the par- 
tial sum, Zj_ 1} with the next value Z^. This use of subscripts brings us 
to our next topic. 



The purpose of subscripting in mathematics is to simplify notation and at 
the same time to remove ambiguity of meaning. In computer work the 
subscript is used in two ways: 

1. Space wise — to specify elements of arrays such as: 

A l A ll A 12- • • • A 1N 

A 2 A 21 A 22 .... A 2N 

A 3 or A 31 A 32' * * * A 3N 



A M A M1 A M2- • • A MN 

This allows reference to elements of an array through simple manipulation 
of the subscripts — just as in normal mathematical thought. 

2. Timewise — to specify the chronological order in which a procedure 
occurs. For example, in Figure 2, diagram B, the subscript i is used 
not only to denote which member of the X and Y array is being operated 
on, but also to indicate how many times the computational step Z = Z + 
(X. - Y.) has been performed. This connotation and use is not so familiar 
to the engineer. The importance of the timewise aspect will become clear 
when iteration is discussed. 
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Absolute Value 



The absolute value becomes important because of the way in which a com- 
puter makes logical decisions . Computer decisions are made on the basis 
of the size of numbers — or, more specifically, upon whether a number is 
positive, zero or negative. For example, in order to do one of two things 
dependent upon whether A is (1) less than 500 or (2) equal to or greater 
than 500, the first thing to do is to subtract 500 from A: 

A - 500 = E 

Then the decision is based upon the size of E. Specifically, if E is neg- 
ative, procedure 1 is done, and if E is positive or zero, procedure 2 is 
done. 

If the difference represents the error in a procedure (that is, 500 is the 
true value and A is the estimate), the absolute value of the error must be 
less than a prescribed amount e and the test is: 

if I E | - e > do procedure 1 

or if |E | - e < do procedure 2 

The absolute value must be used, for in general there is no prior knowl- 
edge of whether the difference A-500 will be positive or negative . All 
the alternatives must be spelled out to the computer in its program. 

In many practical cases the "relative" error is a better measure of the 
error than the absolute error of a result. The relative error test for 
the above situation would be stated: 



if 

if 

Floating Point Arithmetic 



E 



500 
E 



500 



- e > do procedure 1 

- e < do procedure 2 



Engineering calculations often involve handling of very large numbers 
such as Young's modulus, which is 30,000,000 lbs. per sq. in. The 
common procedure is to write such a number as 3.0 x 10? to facilitate 
computation. 

A similar situation exists in handling very large or very small numbers 
in a computer. To get away from carrying a great many digits and to 
eliminate the effort of keeping track of the location of the decimal point 
for these numbers, a floating point notation is used. A common procedure 
is to maintain seven or eight most significant digits of a number plus a 
two-digit "characteristic" to indicate the proper position of the decimal 
point. The characteristic is developed from the exponent of 10 (assuming 
use of the decimal system) . 

For example, the number 

- 00000061957533 
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Errors 



can be represented as 

- .61957533xl0" 6 

However, this representation requires carrying two signs: one for the 
exponent and one for the fraction. This is inconvenient for computers which 
have only one sign associated with a storage location. Therefore, a com- 
mon practice is to add the exponent to some positive base. 

Using a base of 50, the characteristic becomes 50 + (-06) = 44 and the 
internal computer representation of the number becomes: 

4461957533- 

A FORTRAN programmer need not be concerned with this internal rep- 
resentation. The number might appear in a FORTRAN program as 

- .61957533E-6 



There are several sources of errors in computation which deserve con- 
sideration in any problem. 

Initial error: If x is the true value of a data reading and x* is the reading 
used in computation (reflecting an error in measurement, perhaps), the 
initial error is x - x*. 

Rounding error: This type of error results when the less significant digits 
of a quantity are deleted and a rule of correction is applied to the remain- 
ing part. For example, pi, 3. 14159265. . . , rounded to four decimals, 
is 3.1416. 

Truncation error: To simply chop off at four decimal places for pi, giving 
3.1415, would result in a truncation error. Another common source of 
truncation error is in chopping off all terms in an infinite series expansion 
after a particular term. For example, cutting the series for e x at 

e x = 1 + x +x 2 + x 3 
2! 3! 

gives truncation error — sometimes called residual error for series 
approximations . 

Propagated error: If x is the true value of a variable and x* the value 
used in computation, then f(x) - f(x*) is the propagated error. Errors 
may build up in computation ! 



Recursion Formulas 



Very often computation of a set of numbers can be simplified by formulas 
which give one member of the set in terms of other members of the set. 
An example of this is finding the new coefficients of a polynomial when the 
original polynomial has been reduced by the extraction of a root. If r is 
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one root of the polynomial equation a n x n + a -, x 11 " 1 + a^_o x n " 2 + . . . 
a o = 0, then the original 
the reduced polynomial: 



a o = , then the original polynomial can be divided by (x - r) to obtain 



a n xn_1 + < a n-l + ra n> xD ~ 2 + [ a n-2 + r < a n-l + ra n>] x 

z-rJa. x n + a x 11 " 1 + a x 11 " 2 +. . . 
v n n-1 n-2 



n-3 



+. 



a x n - ra x 11 " 1 
n n 



< a n-l + ra n> xn_1 + 1-2 xI1 " 2 + ' ' ' ' 
(a n _ 1 + ra n ) x 11 " 1 - r (a n _ 1 + ra n ) x n " 2 



[ a n-2 + r < a n-l + ra n>] x *~ 2 + 



k k-1 
Let: c^x + c^.jx + + c-jx + c Q = 



be the resulting polynomial (where k = n-1). 
ting coi 

k = a n 



T)'0 T- ll I -rn rs Then, equating coefficients gives: 
-*- " ' s ' "* « I (J CJ g> y | 



CCl) =" A (l-hi) -f 9 M C(l+ l) <>k-l = (a n _! + raj = a^ + rc k 
(^ ^ ' J B c k-2 = [ a n-2 + r < a n-l + ra n >] = a n-2 + rc k-l 



or in general 



c. = a. + re. „ 
l l+l l+l 



Algorithms 



This is a recursion formula for c.. This formula has obvious advantages 
for programming the calculation of the new coefficients . Other examples 
arise in consideration of series expansion, Bessel functions, Legendre 
polynomials and other special functions of mathematics. 



An algorithm is a theorem which may state that a solution to a problem 
exists, and/or a procedure for obtaining the solution. The term is fre- 
quently encountered in computer literature because the form of an equa- 
tion is often all-important in programming efficiently. 

Two examples will be given here: 

1 . To compute the value of the polynomial 

f (x) = a Q x 5 + a 1 x 4 + a g x 3 + a„ x 2 + a. x + a 5 
computation is simplified if it is written as 

f(x) = (((((a ) x + sl ± ) x + a 2 ) x + a 3 ) x + a 4 ) x + a 5 

which indicates that the polynomial may be computed by repeating 
the steps "multiply and add" five times. 
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2. The square root of a positive real number, A, can be computed by 
using the algorithm 

x i+i =i< x i + 4? 

This example also introduces the idea of iteration. What is meant 
by this equation is that to obtain an estimate of A , start with a first 
guess x Q . Substitute it in the right side of the equation to obtain the 
next estimate of A* . A few of the steps are: 

x 2=-K + t> 

x 3=l< x 2 + 4 2 > 



x. , -, = — (x. + ) 

i+l 2 v i x. 
l 

Note that the algorithm states a procedure for solution of the prob- 
lem, not just one formula evaluation. 

Further, note the timewise use of the subscript i; i + 1 indicates a 
result dependent upon the previous result subscripted by i. 

Taking a value of A (say 25) for which the square root is known, and 
performing the indicated operations will make this algorithm clear. 
If 2 is used as a starting estimate, then; 



X-2 



25 ■ "25 



'DO WH!^; ■■. Xa .i-( 2+ f)=7. 



\ ~~ '" ~. ■'••' > l / ' V .i A I 1 / 25 \ _ 

' N ' ■-' ',y, vV./y ;x =—(7.25 + ) =5.35 

C- , f -„ Al y» 2 2 \ 7.25/ 



X 3 = T( 5 - 35+ 5f5)^ 5 ' 01 



There are several ways of deciding the point at which to discontinue 
the iteration: 

a. When differences between successive estimates are less than 
some £ of accuracy | x. +1 -x^ | < £ 

b. When relative differences between successive estimates are 

less than some £ of accuracy 

x; 

i 
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This makes the choice of £ independent of the size of the original 
function value A. 

When relative differences between the function of the estimate and 
the original function value A are less than some £ of accuracy 



(x i+1 ) -A 



This eliminates the possibility of an erroneous solution in a sit- 
uation where convergence is slow and is the preferred type of 
test. 

The general statement for this test in finding a value of x which 
satisfies 



f (x) 



is 



f (x i+1 ) - A 



A 



In succeeding chapters the problem of solving for values of x which 
satisfy equations of the form x = f (x) are discussed. Methods a and 
b remain the same but the test equation corresponding to c above is 



x i+i " f (x,^) 



1+1 



i+r 



t < Vl > 



Complex Number Calculations 



Complex numbers can create some questions of handling and yet are of 
great importance to electrical and electronic engineers . The following 
rules handle complex arithmetic: 

A. Addition 

(a+bi) + (c+di) = (a+c) + (b+d) i 

B. Multiplication 

(a+bi) (c+di) = (ac-bd) + (ad+bc) i 

C . Division 

a+bi = (a+bi) (c-di) = (ac+bd) + (bc-ad) j 
c+di (c+di) (c-di) (c 2 +d 2 ) (c 2 +d 2 ) 
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D. Square Root 

If we let u+vi = (a+bi) 2 

squaring both sides gives u 2 + 2uvi - v 2 = a+bi 

Equating the real and imaginary parts 

(1) u 2 - v 2 = a 

(2) 2uv = b or v = J^- 

x ' 2u 



then substituting from (2) into (1) 

bi. 
4u2 



u 2 - b 2 - a = 



o 
and multiplying through by u gives 

u 4 - au 2 - b 2 = 

T 

and solving the quadratic gives 
■ 2_„ ,r, n >2, /M 2ii 



'-^[(STI*)] 

Only the positive sign is considered since u is real; thus 



and 



wu -[* + ((C + (*) 

(4)v =[r(&) 2 if) 



2U- * 



2\i^ 



If a is positive, compute u by (3) and use (2) to find v. 

If a is negative, compute v by (4) and use (2) to find u. 

E. For a general power of a complex number use 

(a+bi) n = rV 11 " 6 " where r = (a 2 +b 2 ) T 

^=tan _1 A 
a 

and e in = (cos 6- + i sin 6-) = cos n-&+ i sin n-0- 

then (a+bi) n = r n (cos n-0- + i sin n-9- ) 

Transformations and Identities 

Recognition of mathematical transformations and identities which are 
possible may simplify computation. 
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For example, assume in a particular problem that the value for the tan- 
gent of an angle is given: 

tan -9- = A 

and that cos -9- is to be computed. One approach would be to find 
•©• = tan"! A and then compute the cos -G- . However, the cosine may 
immediately be computed as 



cos -9- = 1 



nITTa^" 




Exercises 



saving the evaluation of two trigonometric functions . 

A useful transformation in computing is the logarithmic transformation. 
Suppose we have an exponential function y = k x a . If we take the natural 
log of both sides, we have 

log y = a log x + log k 

thus reducing an exponential relationship between x and y to a linear rela- 
tionship between the transformed variables log x and log y. 

Most of the mathematical techniques which have been outlined here are 
quite familiar to the average engineer. However, they are given this 
emphasis because recognition of their use in particular situations can 
greatly simplify solution of computer problems . 



1. For the set of tabular data shown, construct the difference table to 
determine to what degree of x is y dependent. 









1 


3 


2 


18 


3 


57 


4 


132 


5 


255 


6 


438 



2. Divide 8 + 12i by 3 + 2i. 

3. Using the iteration formula for finding the square root of a number, 
find the square root of 12 accurate to 7 places (the absolute difference 
between two successive estimates of >Il2 is less than .00000005). 



*4. 
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Using only the functions SIN F(X), C(J)SF(X), ATANF(X), EXPF(X), 
L(J)GF(X) and SQRTF(X), write the FORTRAN statement for computing 
all trigonometric, inverse trigonometric and hyperbolic functions. 



5. Construct a block diagram chart to show the logic of computing the 
value of the polynomial 

IT- pq s y' f (x) = a 5 x 5 + a^x 4 + a 3 x 3 + a^x 2 + a-jx + a Q 

where the values of the coefficients a.. , a , a are to be 

read into the computer and the value of f (x) is to be punched out. 

*6. Repeat 5 and generalize the block diagram to show the logic of han- 
dling any degree polynomial. That is, block-diagram the computation 
of f (x) , where 

f (x) = a.]X n + a 2 x n + a n+1 

7. Write a FORTRAN program for exercise 6. 

8. The block diagram in Figure 2 is incomplete. Why? 
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SECTION IV: PRINCIPLES OF ITERATION 



Arithmetic, Mathematics and the Computer 

Although arithmetic is one class of mathematical operation, for purposes 
of clarity it is best to differentiate arithmetic operations from the other 
mathematical operations in discussing computers. Arithmetic is performed 
when any of the operations of addition, subtraction, multiplication and 
division are applied to a pair of numbers . Other mathematical operations 
are performed when one of a variety of operations such as algebraic sub- 
stitution, integration, differentiation, expansion to series representation, 
etc. , is applied to one or more symbols. In essence arithmetic deals 
with numbers while mathematics is the manipulation of symbols. 

Computers are designed to efficiently perform all the arithmetic opera- 
tions . Most computers can also handle in a limited way the alphabetic 
and some special characters which might be used as symbols (as is obvious 
from the use of the computer to translate FORTRAN programs). 

The computer, as designed, cannot do original mathematics. However, 
once the mathematics of a problem are accomplished, the computer is 
most useful in obtaining the required numerical values by performing the 
necessary arithmetic operations. 

In addition, once a particular mathematical problem has been properly 
solved and the program for arithmetic evaluation on a computer has been 
completed, all similar problems, differing only in data values, can be 
henceforth evaluated by the one program . Two examples will illustrate 
these facts: 

A. Given the two simultaneous equations 

a ll x l + a 12 x 2 = b l (1) 

a 21 X l + a 22 X 2 = b 2 (2) 

find values of x and x which satisfy these equations when 

a il =2, a i2 =2 ' a 21 =3 ' a 22 =4, b l =5 > b 2 =6 
To solve this problem we have distinct mathematical and arithmetic steps . 

Mathematically: 

Algebraically using (1) write x.. as a function of x 2 

D-j "" a-« nXn 

x l = f l < x 2> = a & 



11 

Substitute this function f-, (x 2 ) in the second equation to obtain an equa- 
tion in x 2 
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" a 21 a 12 x 2 + a 21 b l + a 22 x 2 = b 2 



(4) 



11 



11 



This equation may then be treated algebraically to give x~ as a function 
of the coefficients only 



)/h 



x 2 = [ b 2 - a 21 b x \ I [ a 22 - a 2i a i2 



(5) 



Y(z)- fe/'O- A(VK 



li 



ii 



Arithmetically: 



Substitute the values given for the coefficients in equation (5) to deter- 
mine x„ 



x = 6- 3. 
2 



4-3.2^ =-1.5 
2 



Use this result in equation (3) to find a value for x 



x _ 5-2(-1.5) _ 1 
1 2 



The computer solution to this problem considers only the arithmetic op- 
erations implied by equations (3) and (5) . A FORTRAN program to eval- 
uate these two equations is shown in Figure 3. 
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Figure 3. 

NOTE the following items of interest concerning this program: 

1. Only the solution equations have been programmed. The original 
equations are not of interest to the computer since the mathematics 
must be performed before using the computer. 

2. Equations (3) and (5) are programmed in their symbolic form and actual 
values for the coefficients are read into the computer before compu- 
tation takes place. By programming in this way the computer can now 
evaluate the solution of any set of two linear equations in two variables , 
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in addition to the particular example discussed. In a sense, by our 
programming, we give the computer the capability of doing nonoriginal 
mathematics. 

B. Solve the following differential equation with stated initial conditions 
for y when x is equal to 1: 



dy 

dx- =x y 



7i7 = x y a) 



^o x ^ 






G/VD 



y = 1 when x = 
Mathematically: 

Integrate (1) for y as a function of x 
dy 

dx- = *y 

iS = xdx 

x 2 
lny =y-+k 

y = e * 
Considering the initial condition y = 1 when x = 0, then 
k = 



Arithmetically: 



Since the analytic solution is y = e 2 it is now necessary only to sub- 
stitute the value of x for which y is to be determined and perform the 
indicated arithmetic. 
j_ 
y =e 2 = 1.64872 

A FORTRAN program to accomplish the required calculations is shown 
in Figure 4. 
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Figure 4. 
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NOTE the following: 

1. Again only the analytic solution is programmed. 

2. The program is written so as to allow many values of y to be com- 
puted over any range in x. 



Numerical Analysis 



The two problems chosen happen to be of types for which the mathematical 
analysis is well defined. The variables in question may be equated to 
functions of the known quantities — analogous to the formula y = f (x, z) as 
discussed in the first chapter. 

Various analytic techniques successfully handle a large class of the math- 
ematical problems which arise in the description of the physical world. 
However, a vastly larger number of mathematical problems are not sub- 
ject to the normal analytic procedures of algebra, analytic geometry, or 
the calculus. 

For the theoretical mathematician or physicist it may be sufficient, when 
meeting these problems, to be able to show that one or more approxima- 
tions will allow an analytic solution of the problem to any desired degree 
of accuracy. 

The engineer, on the other hand, must have a solution from which accurate 
numerical results can be obtained if the theory is to be of value. These 
results must be obtained at a reasonable expense of time and effort. A 
mathematician may be content to know that a series will converge . An 
engineer must be concerned with the number of terms in the series re- 
quired for a given accuracy. 

Numerical analysis provides techniques for finding concrete numerical 
results for mathematical problems of the type under discussion. This 
manual describes many physical problems which require numerical analy- 
sis. Some of the simpler analysis techniques are discussed. The reader 
is referred to the many excellent texts available on numerical analysis for 
an adequate discussion of particular analysis theory. 



Algebraic Equation 



Problem A belongs to a very special class of algebraic problems, a set 
of n nonhomogeneous linear equations in n variables. Fortunately, a 
large class of engineering problems such as the force and moment prob- 
lems of stress analysis give rise to this simply handled type of equation. 
The diagram below shows a simple force equation example giving rise to 
a set of two linear equations . 
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IF X = 
lF y = 



T 2 T l 








W 



T, cos6-- T cos (J) = 
1 2 

T 1 sin-0- + T 2 sin (J) = W 

T, and T are readily solved for specified values of -9- , (J) and W. 

Consider a set of equations similar to A but containing a simple 
nonlinearity: 



a ll x l + a 12 x 2 ~ b l 

3,nnXn 

a_ x. +e z ^ z = b. 



(1) 
(2) 



21 1 " "2 

The previous algebraic steps would give 



a 12 ^ bi 

X- = Xo + _±_ 

1 a ll * a n 



(3) 



and: 



a M b« -a 2 2X 2)/ / r . a 21 a 12) 



x 2= <V IJjL" e 
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(4) 



or with appropriate substitutions for (4) 

bxo 
x = ae * + c (5) 

This gives an equation in x 2 only, but from the equation it is impossible 
to directly evaluate x 2 . One approach is to expand e bx 2 to several terms 
of its series: 

e bx 2 = i + b^2 + (^ + (b^2) 3 + . . . . 
1 2! 3! 

However, for sufficient accuracy it may be necessary to carry a large 
number of terms in the series resulting in a polynomial equation being of 
degree greater than 4. In such a case, the equation is again impossible 
to evaluate directly. This chapter will discuss handling equation (5) and 
polynomials by means of iteration. 
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Equations of Calculus 



Similarly the differential equation of example B has an analytic solution. 
However, the differential equation 



dy 
dx 



x 2 + y 2 



has no corresponding analytic solution. To find values of y for corre- 
sponding values of x, a numerical approach is necessary. Likewise, 
many integrals may be obtained only by numerical approximations. Again, 
some simple techniques will be given for handling this type of problem. 



Iteration Type Problems 



There are three principal reasons for the use of iteration with computer 
handling of engineering problems. 

1. The mathematical statement of the physical problem requires an itera- 
tion approach for evaluating one or more of the variables. 



The chart below shows examples of four types of mathematical problem statements. 



Single Equations 


Multiple Equations 




JL x/2 
x= 5 e 


M H =A H- m 
M He = A He " 2m 
_ R He- R H 


Direct Iteration 
Sufficient 




R H/M H - R He /MHe 




(J) =tan(J)-K 


B = 90 - o - sin 1 J 

L = A sin (90° - B - TS -<<) - (R-E) 


Numerical Techniques 
Required 



2. Many designs are to be evaluated in a search for the best design, that 
is, optimization of design. Often this problem can be reduced to the 
previous problem by selection of an appropriately expressed mathemat- 
ical criterion of what is optimal. 

3. The mathematical expression for the physical problem is to be eval- 
uated for many values of one or more known parameters — such as 
time, in a problem of motion, or degree of rotation, in a geometry 
problem . 



A. One problem which illustrates iteration is the solution of an equation 
which arises frequently in absorption problems of optics, electric 
fields and nuclear engineering: 



bx 



(1) 
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where a and b are constants. 



"DO loi-MLe, 



This cannot be solved explicitly for x, so the following iteration pro- 
cedure is used: 

1. Make a guess at x. 

2. Use this in right-hand side of equation (1) to give a new value for x. 

3. Call this new value of x the next guess. 

4. Repeat steps 2 and 3 until two successive guesses either agree or 
differ by an amount less than the allowable error. 

A FORTRAN program to solve such a problem, where 

x = 0.2e 0,5x 

is as follows: 



READ 10, A, B, X 

1 PRINT 20, X 

2 XNEW=A*EXPF(B*X) 

TEST=ABSF (X-XNE W) 
~"TF(TEST-700006)4 T 4A 

3 X=XNEW 
•GO-TQ-L- 

4 PRINT 20,XNEW 
END 






Read values for A(0.2), 
B(0.5) and initial estimate 
of x (1.0). 

Print estimate of x. 

FORTRAN arithmetic state- 
ment of equation to be solved. 

Find absolute value of dif- 
ference in last two estimates. 

If difference is less than or 
equal to .00005, go to step 
4. Otherwise, go to step 3. 

Store new estimate of x. 

Return to step 1 and repeat. 

Print last estimate of x. 

End of program . 



If this were translated into machine language and the resulting program 
run on a computer, the successive estimates of x would be: 

X 

1.000000 

.329744 

.235848 

.225032 

.223818 

.223682 

.223667 
40 



The last value in the list satisfies, the requirement that the difference 
between it and the previous estimate be less than 0.00005 in absolute 
value. Hence it is the solution. 

B. The simple iteration procedure in problem A seems like a good method 
to solve transcendental equations — and in many cases it is. However, 
on application to the following problem — much like problem A except 
that a constant has been changed — it fails: 

f(x) = x - Ae Bx = A = .2, B=3 (1) 

The reason for the failure of convergence of the direct iteration technique 
should be explained. Direct iteration has the formula 

x _= f (x ) 
n+1 n' 

This indicated procedure will converge only for 

-Kf'(x )< 1 
n 

Bxn 
where f (x ) = A e for this problem, and the prime indicates the first 

derivative. The choice of a larger B makes f (x) > 1 and therefore makes 

convergence by the direct procedure impossible. 

Since simple iteration fails to produce convergence , a different method to 
obtain an estimate must be found . The Newton-Raphson method answers 
this need. The Newton-Raphson iterative equation is: 

F(x n ) 

x n+l ~ x n~ 

F'(x n ) 

Applying this iteration scheme to equation (1): 

Bx n 
F (x^ = x n - Ae 

and 

Bx n 
F' (Xjj) = 1 -BAe 

Applying the Newton-Raphson formula: 

Bx n 
x n - Ae 



x n+l = *n " Bx n 

1 - BAe 

To solve this problem in FORTRAN, replace statement 2 in the previous 
FORTRAN program with the arithmetic statement for (3): 

2 XNEW=X- (X-A*EXPF(B*X)) / (1.0-A*B*EXPF(B*X)) 

If this FORTRAN program were translated into a machine-language pro- 
gram and run on a computer, the successive values of x would be: 
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1.00000 
.19742 
.22364 
.22366 

Note the increased rate of convergence obtained using this new method; 
only three new estimates are computed, compared with seven in the pre- 
vious example. 

C . The following engineering problem will illustrate one origin of the 
type of equation just discussed. The analysis makes use of the clas- 
sical optimization principle of equating the derivative of a function of 
one variable to zero. The values of the variable which satisfy the 
resulting equation are those for which the original function is either 
a maximum or a minimum. (Since engineers are often interested in 
optimization of design, it is unfortunate that so few design problems 
can be expressed in a form suitable for the use of this' principle . ) 

For high-potential conduction through walls, tubular insulators are used 
which are covered with metal on the inner and outer surfaces. 



Bore 




» r 



What must be the ratio of the external diameter, 2R, to the bore width, 
2r, for the cross section Q to be a minimum? 

The ratio q of the line voltage to the maximum admissible field strength 
is given by 



q - r ln(R/r) 
The cross section Q is given by 

9 9 

Q = 7T (IT - r ) 
According to (1) , letting x = R/r , then 
r = q/ln x 



(1) 



(2) 



or 



Q=7rq 2 (x 2 -l)/(lnx) 2 



(3) 
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Since Q is to be a minimum, it is necessary to find a point on the curve 
defined by (3) where the slope is zero. To do this the first derivative is 
set equal to zero: 



dQ 
dx 



= it Q 



2x 



2(x z - 1 ) 
(lnx) 2 x(lnx) 3 







(4) 



From the bracketed expression: 



x - (x - l)/x In x = 



or 



In x = 1 - \/yl 

This in turn gives 

1-1/x 2 

x - e ' = 



(6) 



The solution for x is readily recognizable as the type of problem just 
discussed. 

An approximate graph of this equation will illustrate some interesting 
points: 



.5 








/(x) 










.5 


l.0\ 


1.5 J^ 2 


-.5 









An examination of the equation (6) will show that, as expected from the 
graph, an initial estimate of 1.0 will give back a solution value of 1.0 using 
either the direct iterative procedure or the Newton-Raphson technique . 
However, from the problem definition this is a trivial solution (equal inner 
and outer radii will give no cross section) and the root of the equation 
desired is the second one at 2.2. 

For ah initial estimate chosen beyond the minimum point, indicated by 
"a", the techniques will converge to the proper solution. An explanation 
for this will be given in the next problem to be discussed. 

D. Another example can be taken from an equation which occurs when 
working with helical gears. This is 



K = involute of = tan (J) - d) 



(1) 



where K is found in a table and it is necessary to find an angle (J) satis- 
fying the equation. 

The Newton-Raphson method can again be used but it is worthwhile to stop 
and consider a good way to get a first estimate of (J) . Since the rate of 
convergence — and in some cases, the possibility of converging — will 
depend on the initial guess, this is a subject worthy of some thought. 
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The following graphic presentation of this problem will clarify the situa- 
tion. 

Plotting Newton's function: F(0) = tan <j) - (j) - K gives the curve as shown: 




Possible range of 
• initial estimate 
for convergence 



Enlarging the area of intersection of F ((})) with the <j) axis shows the basis 
for Newton's procedure. 




If (J) is the first estimate of the root, the figure indicates that the next 
good estimate would be the intersection of the projected slope of F ( (f> ) 
evaluated for (J) n . The following use of the construction angle -6- shows 
this to be precisely Newton's procedure: 



or 



tan -9- = 


= F((J) 
v n 


n) 
■On+1 


n+l 


= <D 

n 


-F«D n ) 
*"(<P n ) 



= F' (<J) n ) 
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The dotted lines in the figure show the limits on the range of the initial 
estimate for which Newton's technique will converge to the proper value. 

A good first estimate can often be determined by a careful examination of 
the problem to be solved. In this case the first two terms of the trigono- 
metric series for tan (J) give: 

~. <D 3 
tan = (J) + — (2) 

Substituting this in (1) and solving for (}) gives 

(J) = (3K) 1/3 

Thus an initial guess which should be reasonably close to the solution can 
be computed from the given value of K. This can be an important advan- 
tage where convergence is likely to be slow. 

To illustrate the effect, the problem has been solved for two values of 
K (0.001 and 0.01). In each case, the initial estimate of (J) = 0.52359 
radians (30°) was used rather than a value computed as described above. 
The results obtained are as follows: 

<t)(K=.001) 0(K=.01) 



.52359 .52359 

.36535 .39235 

.25482 .32546 

.18624 .30818 

.15296 .30702 

.14440 .30701 
. 14385 
. 14372 

Note that for K = .01 fewer iterations were required because the initial 
estimate was closer to the correct result. Had the approximation formula 
developed above been used, the initial estimates would have been 0. 14423 
and 0.31072 respectively, and no more than two iterations would have been 
required in either case. 

E. A simple geometry problem which lends itself to iterative solution is 
that of finding the angle for which the arc and chord in a circle of 
given radius will enclose a stated area. 
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Area A is given by 



2 o 

A = tt r ^6- _ r sin 6- 

360 2 



(1) 



Since -0- cannot be solved for directly, the Newton-Raphson technique is 
used to find a solution. Converting -9- to radians and applying the N-R 
formula gives the iteration equation: 



where 



■0 n =-& _ F(-0- ) 
D+1 " TTT^J 

F(*) = r2 *„ " ^ sin * n " A 



and 



F'(-e-) =if. - r 2 cos-9- n 



It is interesting to note how frequently an iterative problem is given rise 
to by asking an unfamiliar question to a familiar relationship. Generally 
equation (1) is used to determine the area A. In a sense, iteration arises 
when the "answer" is given and the "question" is to be determined. 

F. Very often iteration can be the most convenient procedure in finding 
a root of a cubic such as 

ax +bx +cx + d = 

Newton's technique will generally suffice here. 

For polynomials of degree greater than 4, some type of iterative proce- 
dure must be used. There are many techniques available to handle the 
higher-degree polynomials, but they will not be discussed here. One of 
the important advantages of the FORTRAN language is that it makes it 
more readily possible for computing organizations to establish standard 
programs for such problems and exchange them. 

G. This problem is an example of how a simple circuit can give rise to a 
cubic equation. Generally it is simpler to make successive approx- 
imations than to solve the cubic . 



Nonlinear 

Passive 

Component 

I-kV 3 ' 2 

k = 10 amps/ volts 3 / 2 



T 



Linear 
Active 
Component 

Emf = 100 volts 

R = 5000 ohms 
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Problem: Given the circuit shown, find the potential difference V from 
(1) to (2) and the current I flowing in the circuit. 

Using the basic relationship for the potential drop, 

V = E - IR 

and considering the active component, gives 

V = 100 - 5000kV 3 ' 2 
or 

V = 100 - .05V 3/2 

a cubic equation easily solved using Newton's technique. 

H. An interesting example of the use of iteration can be taken from the 
computations employed in spectroscopy. Spectroscopy and astronomy 
are two areas of science in which it is meaningful to speak of numbers 
accurate to seven or eight significant figures . The success of iterative 
approximations in the field of spectroscopy demonstrates beyond ques- 
tion its possibilities for accuracy. 

Given experimentally determined values of the Rydberg constants for 
hydrogen (R H ) and helium (Rjj e ) and the atomic mass (A H , A R ) of 
the elements, it is possible to compute the mass of the electron. 

Assume the following data is given 

R R = 109,677.68 cm" 1 A R = 1.00812 

R Re = 109,722.34 cm" 1 A =4.00388 

if R^ is the Rydberg constant assuming the mass of the nucleus to be 
infinite, the following relationships exist where m is the mass of the 
electron: 

R co = ^X> % (^ 

R = (1+™ ) R 

00 M He He (2) 

The atomic mass includes mass of the electron: 



A H . = M R + m (3) 

A TJ = M rr + 2m 

He He 
This gives a useful first approximation 



M H =A R (4) 



M„ = A~ 
He He 
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Equating (1) and (2) gives a first evaluation of m: 
R„ R. 



m 



= He - H (5) 

R H R He 



/•Use 



5?- O 



M M 
H He 

To obtain subsequent approximations, equation (3) is rewritten as: 

M„ = A„ - m (6a) 

and 

M He = A He- 2m < 6b > 

and, using the new values m is re-evaluated by (5). 
• A FORTRAN program for this procedure is: 

20 EMASS-0.0 

30 RH=109677.68 

.'"'■ = loiL-n.gft 40 RHE=109722.34 

RNe=/6<J72'Z.-3i/ 50 AH=1. 00812 

/-) H = J . o oS'j 7 6 AH E=4 . 00388 

iy-!e "w a, D p-7c.c 70 EPSLN=. 00000005 

' a ' ^ c p>l £MH=AH~EMASS 

80 EMHE-AHE-2.0#EMASS 

90 J:MASP=(RHE-RH)/(RH/EMH--RHE/£MHE) 

100 PUNCH 20,EMASP,EMH,EMHE 

K»h i HO IF(ABSF(EMASS-EMASP)-EPSLN)3,3,2 

>-> LOiiiir/ , L ' 2 $MASS=$MASP 

''/^~" U,5s ^ 120 GO TO 1 

M ^ J -M A5 p - S -^ > ). 3 STOP 

^ ) 
Running this program on a computer produces the following results: 

s ^H ^m 

. 00054871426 1. 0081200 4. 0038800 

.00054836568 1.0075713 4.0027826 

.00054836592 1.0075716 4.0027833 

It is significant to note that the simplest iterative technique produces the 
necessary accuracy for this application in only three iterations. 

I. Earlier it was stated that iteration could be used in obtaining numerical 
solutions to differential equations. To illustrate, consider a problem 
for which the exact integral is known: 

and find y as x varies from to 1.0. The initial conditions are that for 
x = 0, y = 1.0. 
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One method of evaluating the equation numerically for y = f (x) over the 
range in x (Euler's method) is illustrated in the figure below: 



n-H 




Numerical 
Approximation 



hAxH 



Taylor's series may be used to estimate the value of a function y = f (x) in 
the vicinity of a given point y 



y n +i 



y * + (4^ + (&)^ 






Consideration of the first two terms only gives 



(i) y n +i 

where 



" '» + Wn ** 



< 2 > 1&4 ■ v. 



from the given differential equation. 
From the figure it is clear that 
(3) x n+1 = x n + Ax 



Successive evaluation of equations (1), (2) and (3) forms the iterative 
scheme necessary to repeatedly calculate the numerical approximation 
for y at regular intervals in x over the required range in x. 
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To get started in the use of these iterative equations, the initial conditions 
are applied as follows: 



$) = **> 



i^uc - 



•6 5 



Do 



w , TX-LX 

X-= O BY >C-5 TO | 



x± = x Q + Ax 



The FORTRAN program to evaluate the stated problem, and a sample calcu- 
lated result follow: 



3 DELX =.05 

4 X=0.0 

5 Y«1.0 

1 PRINT 20,X,Y 

6 DELY=X-ttY-*DELX 

7 X=X+DELX 

8 Y=Y+DELY 

9 IF(X-1.0)1,1,2 

2 STOP 
END 

Taking Ax = . 05, the following solution was obtained 



. 00000 


1.00000 


. 05000 


1.00000 


. 10000 


1.00250 


. 15000 


1.00751 


. 20000 


1.01507 


. 25000 


1.02522 


. 30000 


1.03803 


. 3^5000 


1.05361 


. 40000 


1.07204 


. 45000 


1.09348 


.50000 


1.11809 


. 55000 


1.14604 


. 60000 


1.17756 


. 65000 


1.21288 


. 70000 


1.25230 


. 75000 


1.29613 


. 80000 


1.34474 


. 85000 


1.39853 


.90000 


1.45796 


. 95000 


1.52357 


1.00000 


1.59594 
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The analytic solution to (1) for the stated boundary condition is 

2 

y = e 



and for x = 1.0, y = 1.64872. The error in the numerical solution is 
1.64872 - 1.59594 = .05278, or about 3%. Considering the size of the 
interval, this is not out of line. Increased accuracy could be obtained 
by decreasing Ax or by using a more sophisticated technique — and there 
are many. 

Consideration of the first alternative should take into account that de- 
creasing the value of A x may increase the rounding error and hence make 
it necessary to carry along a greater number of significant digits. The 
second alternative must be viewed in the light of the accuracy required. 
The more sophisticated techniques will require greater programming 
effort and more computer time. It is not reasonable to spend this time 
and effort to obtain accuracy of .01% if the required accuracy is only 2% ■— 
particularly in cases where the input data is accurate to, say, 5%. These 
points should be considered in programming any problem for a computer 
if the most economical use is to be realized. 

The more sophisticated techniques usually make use of a weighting of 
several previous estimates for the derivative in an equation similar to: 

y n +i = y n + w i < w 2 y'n + w 2y n +i + — > Ax 

where the prime indicates a derivative and the w^ are selected weighting 
constants. 

However, the logic of predicting the next value for y in terms of the last 
calculated point and the derivative evaluated according to the given differ- 
ential equation remains the same. 

J. It is easy to extend the simple numerical integration procedure just 
discussed to higher-order differential equations. For example, sup- 
pose the following is to be solved over the interval to 1.0: 

£_y + a -L + Bxy = 
dx^ dx 

where x = 0, y = 1.0, (dy/dx) = 1, A = 2 and B = 3. 

The analytic solution must be in terms of an infinite series. Therefore, 
whether a series is determined or numerical integration is performed, 
the solution involves considerable computation. The series representa- 
tion allows control of the error while the stepwise integration procedure 
to be described may not. 
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An approach to this problem is to start with 

(S)o =-A«>y/dx> -BV 

and then obtain successive values of each of the variables from 

y n+l = yn + < d y /dx >n Ax 

Vl =x n +Ax 

(dy/dx) Q+1 = (dy/dx) n + (d 2 y/dx 2 ) n Ax 

(d 2 y/4x 2 ) n+1 = -A(dy/dx) n+1 -Bx n+1 y n+1 

A computer solution to this problem is shown in Figure 5. 

* * * * 

In the first few problems just presented, iteration was used to improve 
the estimate of a single solution. In the case of the differential equations, 
iteration provided successive values of y for finite changes in x. These 
are two entirely different concepts of iteration. 

1 A-2.0 

2 B=3.0 

3 XN s 0.0 

4 YN=1.0 

5 DYDX=1.0 

6 D2YDX S -A*DYDX-B*XN#YN 

7 PRINT 20, XN,YN,DYDX,D2YDX 

8 IF(XN-1.0>9,14,14 

9 YN=YN+DYDX#.05 

10 XN^XN+.OS 

11 DYDX-DYDX+D2YDX* .05 

12 D2YDX=-A*DYDX-B*XN*YN 

13 GOTO 7 

14 STOP 
END 



XN 



YN 



DYDX 



D2YDX 



00000000 


1. 


.0000000 


1.00000000 


-2, 


.0000000 


05000000 


1, 


.0500000 


.90000000 


-1, 


.9575000 


10000000 


1, 


.0950000 


.80212500 


-1, 


.9327500 


15000000 


1, 


.1351063 


.70548750 


-1, 


.9217729 


20000000 


1, 


.1703807 


.60939886 


-1, 


.9210261 


25000000 


1, 


.2008506 


.51334756 


-1 


.9273331 


30000000 


1, 


.2265180 


.41698091 


-1 


.9378280 


35000000 


1, 


.2473670 


.32008951 


-1 


.9499144 


40000000 


1, 


.2633715 


.22259379 


-1 


.9612334 


45000000 


1, 


.2745012 


.12453212 


-1 


.9696408 


50000000 


1, 


.2807278 


.02605008 


-1 


.9731919 



Figure 5 Continued on top of next page 
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XN 



YN 



DYDX 



D2YDX 



Exercises 



.55000000 


1 


.2820303 


-.07260951 


-1 


.9701310 


.60000000 


1 


.2783998 


-.17111607 


~l 


.9588875 


.65000000 


1 


.2698440 


-.26906045 


-1 


.9380749 


.70000000 


1 


.2563910 


-.36596420 


-1 


.9064927 


.75000000 


1 


.2380928 


-.46128884 


-1 


.8631311 


.80000000 


1 


.2150284 


-.55444540 


-1 


.8071774 


.85000000 


1 


.1873061 


-.64480427 


-1 


.7380221 


.90000000 


1 


.1550659 


-.73170538 


-1 


.6552671 


.95000000 


1 


.1184806 


-.81446874 


-1 


.5587323 


1.00000000 


1 


.0777572 


-.89240536 


-1 


.4484609 



Figure 5 

Engineering computational problems may be classified in the following 
way — the first category being of trivial interest here, although certainly 
not trivial as far as computing is concerned: 

1. Parameters specified, explicit equations. 

2. Parameters specified, iteration required. 

3. Parameters variable, explicit equations.. 

4. Parameters variable, iteration required. 

Problems A through H are typical of the second category. Differential 
equations may fit in either the third or fourth category — examples I and 
J falling in category 3. The next subject to be discussed — mathematical 
models — also comes under either category 3 or 4. 



1. For problem C, find Newton's iterative equation and write the 
FORTRAN program to solve the problem. 

2. Write Newton's equation to find a root of 

5x 3 + 2x 2 + 3x + 1 = 

3. The following is a problem which occurs in the planetary gear systems 
currently used in automatic transmissions on automobiles. The curves 
y = sin cut and y = e~ at are plotted on the same set of y, t axes. De- 
termine the smallest positive value of t at whteh they intersect. 
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a. Set up the equation for which a root is to be determined. 

b. Write Newton's iterative equation to find the root. 

*4. A common numerical procedure for calculating integrals is the use of 
Simpson's rule which may be stated: Given y as some function of x, 
the integral of the function y over the range x = x , to x = x may be 
approximated by: n 



x =x. 



n y dx s^L. (yi + 4y2 + 2 y3 + 4 y4 + 2yn _ 2 + 4^.! + y n ) 



\jx = x^ 
where x^ = x^_- L +Ax 

and y. is evaluated for the corresponding x. . 
Write a FORTRAN-language program to evaluate 
x=1.0 1 



S 



y dx where y = 1+x 2 

x = 



using a Ax = .05. 

5. The "spectroscopy" problem can be solved directly. Perform the 
analysis and write the FORTRAN program to accomplish it. 
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SECTION V: MATHEMATICAL MODELS 



The suggestion was offered in the introduction that the "build and try" ap- 
proach to engineering is being replaced by "construct mathematical model 
and simulate". The problems mentioned thus far for computer handling 
certainly would not warrant a "build and try" approach. 

In point of fact, engineers are not presented mathematical problems of this 
type directly, although these problems readily occur as their work develops. 
They are presented with the task of designing or creating a part or system 
which will perform a given function. 

The actual size and complexity of the system may vary from a gear train 
assigned to one engineer, up to an automobile which requires the combined 
design efforts of many engineers. However, the emphasis remains on 
function. 

The engineer is measured in terms of the performance of the design result 
and the time required to accomplish to design. 

The problems considered in this section will illustrate the use of mathe- 
matical models in evaluating the performance of proposed designs. This 
leads to selection of better designs and a reduction in the design time. 

The following chart illustrates the analogy between the "build and try" ap- 
proach and the computer approach. 





Specify 
Function 
and 
Characteristics 




Specify Input 
Values and 
Criteria for 
"Goodness" 






\ 


' 


1 


1 






Build 
Physical 
Models 


Construct 
Mathematical 

and 
Logical Model 










1 


\ 


- 










1 


f 












Test 
Models 


Test 
Models 








Modify 






Modify 






























* 








Build 
Prototype 






1 


' 






Make Final 
Adjustments 
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A few comments concerning this parallel flow of the two approaches are in 
order. 

An outstanding advantage of the computer approach is the construction of 
a prototype as the first physical build. This is a common result of the 
computer approach and has often eliminated the cost of extensive, time- 
consuming physical construction. 

There exist other significant differences. The specification of parameters 
(size, shape, material) may often be presented in qualitative terms for the 
"build and try". These qualitative values must be replaced by specific 
quantitative values for use with a computer. Criteria for a "better" design 
must be expressed in precise mathematical terms. 

The model is much more than mathematics. Mathematics may suffice to 
describe completely a component of the system, but rarely will known 
relationships be adequate to describe completely the intereffects of com- 
ponents working together. Considerable logic, often from engineering ex- 
perience, must be used to create a realistic, accurate representation of 
the system to be studied. 

Finally, in testing the model, possible modifications may be part of the 
mathematical model and programmed so that automatic modifications take 
place as a result of decisions based on steps in the calculation. Of course, 
these modifications must have been considered during construction of the 
model and often represent a searching for a "better" design according to 
the defined criteria. It can be expected that the unforeseen type of basic 
modification to the model will also be required, before the testing is 
completed . 

The problems which follow are selected to illustrate various aspects of 
model building for computer analysis . Most of the descriptions are con- 
cerned only with the facet of the physical problem presented which requires 
something more than conventional analytic techniques. In actual use, many 
more calculations are performed than indicated here. 
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Solenoid Design 



This problem is most typical of design problems where several variables 
are involved and there is no obvious procedure for arriving at the best 
design (category 4 — parameters variable, iteration required). 

Problem Statement: 

Choose a wire for a solenoid such that the power consumed will be less 
than some fixed amount and the solenoid will give a fixed pull. 



K. 



Engineering handbooks can supply the following set of relationships: 



Pull 



IN 



F =- 



(1) where N=total number of turns 



Current 



-I 



(2) where V=voltage applied 



Power 



P = 



(3) 



Resistance 



Wire Length 



R = 4oC 7L 
IT 

L =TT DN 



(4) where/ 7 = resistivity of wire 



(5) 



Solenoid Diameter D = d+t 



(6) 



Coil Thickness t =^- 



N 
Number of Layers n =7jj[~ 



(7) where o< reciprocal of wire 

diameter (wire size) 



(8) 



This is often the way a problem is presented to an engineer. If use of a 
computer is contemplated, the usual rules of analysis apply. That is, it 
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is best to~do-a little algebra before plugging in numbers. Too often the 
tendency is to choose wire size<< , start at equation 8 and work back, 
equation by equation, to arrive at power finally by equation 3. Applica- 
tion of algebra, however, gives: 

FV 
P= (^-« 2d ) W 



4F/ ? je 



or power directly in terms of wire size. Savings in computer time can be 
just as important as savings in engineering analysis time. So it always 
pays to do a little analysis before muddying up the water with numbers . 

The important aspect of this problem is not that one can now compute the 
power output for a given size and thus choose an acceptable wire size, 
but rather that a computer program can be generalized such that, given 
any requirements for a solenoid, many different combinations may be 
tried quickly to find the suitable combination for the job at hand. Thus, 
to design a solenoid for a given pull F, it might be desirable to vary SL , d, 
and V or, perhaps, to limit t. Using equations 1 to 8, a variety of opti- 
mizing programs can be written for solenoid design. The question which 
would determine the general form of each program is: What particular 
specifications are to be met? 



Heat Transfer Problem 



An insulating wall is made of three parallel layers of different insulating 
materials. The outside temperatures, t and t , are known. The con- 
ductivity k. of each layer, a straight line function of the mean temperature 
t. , is known. The problem is to find the quantity of heat Q passing through 



a unit area of the wall per unit of time. 



ti 



*2 



*3 



t4 



x i 



x 2 



« x 3 ► 



Q is given by 



Q = 



ti -t, 



A.-. An 



(1) 
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However, the k. are dependent on the t. in the following manner: 



k-, = a, t. + b 1 



h" a 2^ +h 2 < 2 > 



k 3 = Vs + b 3 



where the a. and bj are given constants for the particular intervals. At 
the same time, the t- are given by 



Vt-Oi-^i^) 

1 1 2~ 2 

t 2 = t 2 -(VS)=i(' 2 +t 3 ) 

2 2 



(3) 



but t ? and t~ are not known. 

If the k. are known, t and t may be determined from the steady-state 
requirement that the quantity of heat passing through each layer per unit 
of time must be the same . 

Or 

k l (h ~ fc 2) = k 2 < fc 2 " *3) = k 3 ^3 " fc 4) 

Xo 





x l 






x 2 




which 


gives 
















x 3 


(h- 


t 4 ) 




'3 = 


V 


x l 
k l 


x 2 


x 3 

+ k 7 




V 


■h- 


x l 
k l 


(h- 


t 4 > 




x l 

k l 


+ f2 

k 2 


+ ^3 
k 3 



(4) 



But this produces a vicious circle; the k. are needed to determine t2 and 
to and vice versa. Thus, iteration is required. 

The procedure for this is as follows: 

1. Make a reasonable guess at t£ and t„. 

2. Use (3) to determine the t^. 

3. Use (2) to determine the kj. 
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4. Use (1) to determine Q. 

5. Since values are known for kj, recompute t 2 and t 3 using equations (4) 

6. Repeat steps 2 through 5 until two successive values for Q agree. 

This problem converges to a solution nicely using the simplest iterative 
scheme, x = f(x,y,a). Experience has shown that the initial estimates of 
the internal temperatures are in no way critical except that t^ must be 
greater than t 2 , which must be greater than t3, which must be greater 
than t4. 

The data used in this sample computation is as follows: 

b. a, 

ing Brick .0623 .00010 

Kaolin Insulating Firebrick . 0255 . 00005 

2.4395 .00060 



Layer 1 


Kaolin Ins 


Layer 2 


Kaolin Ins 


Layer 3 


Magnesite 


t ± = 1400° 




t 4 = 200° 




x = 2.15" 




x 2 = 2.0" 




x 3 = 6.0" 





The first set of results in Figure 6 shows the iterated solution to this 
problem giving the final value for the quantity of heat transmitted in a 
unit_time as_Q = 25.38 Btu/min with the corresponding mean temperature 
tj_, t 2 , and t 3 . An absurd case is also shown by the second set of results 
in Figure 6, with the input data changed in order to illustrate that con- 
vergence is not always possible. The data changed is 

H 
.02 

.01 

-.05 

Actually it can be seen that the reason for the divergence of the second 
set was that the large negative slope a.3 caused k 3 to become negative (an 
impossible situation physically, since the coefficients of conductivity may 
not be negative). The value of -.05 for a 3 is not a reasonable value. It 
was chosen here to indicate that a mathematical procedure of this type is 
valid only within some range of selection of data. 
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PQ=0.0 

1 READ 10, Bl ,B2,B3,T1 ,T4,X1,X2,X3 

2 READ 10, Al ,A2,A3 

3 T2=2.-*(Tl-T4)/3. + T4 

4 T3=(Tl-T4)/3.+T4 

5 T1B=.5*(T1+T2) 

6 T2B=.5-(T2+T3) 

7 T3B=.5-"-(T3+T4) 

8 Z1=A1-*T1B+B1 

9 Z2=A2-"-T2B+B2 
19 Z3=A3*T3B+B3 

11 Q-(T1-T4)/CX1/Z1+X2/Z2+X3/Z3) 

12 PRINT 20,Q,T1B,T2B,T3B 

13 IF(ABSF(PQ-Q>-. 00005)2, 2, 14 

14 PQ = Q 

15 T3=T4+X3/Z3*Q 

16 T2 = Tl-Xl/Zl-"-Q 

17 GO TO 5 

18 END 

Iteration Q "L I2 "£3 

1 26.9~25558 1200.0000 800700000 4007o"0000 

2 25.146838 1241.2234 671.36960 230.14618 

3 25.385972 1254.9914 684.25925 229.26788 

4 25.384732 1254.6856 684.23785 229.55224 

5 25.384376 1254.6691 684.21790 229.54884 

6 25.384383 1254.6698 684.21825 229.54845 

Iteration Q^ T^ T2 "£3 

1 -384754.98 1200.0000 800.00000 400.00000 

2 188995.29 18589.197 83719.950 65930.755 

3 145.56898 853.61840 81.496675 27.878275 

4 38318.358 1390.8673 1208.5344 617.66710 

5 -823.77876 -77.502100 -4718.9752 -3841.4731 

6 4076.2550 804.76095 192.05571 187.29476 

Figure 6. 
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Rocker Arm Cam Problem 



The accompanying figure shows a rocker arm and cam to be used in a 
fuel pump. 




Specified are: 

1. Axis of rotation of the off-center circular cam (X£, Y<£. 

2. Eccentricity, E, and radius, R, of cam. 

3 . Axis of rotation of rocker arm (X, , Y ) . 

4. Required angular rotation, o(. , of rocker arm to give desired travel 
of pull rod. 

To be determined are: 

1. Necessary angle, B, between rocker arm and vertical. 

2. Necessary offset, L, for rocker arm. 

To obtain a solution, "boundary conditions" are first considered to obtain 
expressions for B and L. For the "high" position of the rocker arm 



B = 90° - tf - arcsin ( L+R+E ) 

A 

For the "low" position of the rocker arm 



(1) 
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L = A sin (90° -B- tf -<£ )-(R-E) (2) 

The direct iterative procedure for solution is: 

1. Make a reasonable guess at L. 

2. Determine B from equation (1). 

3. Then use this value of B in (2) to obtain a value for L. 

4. Repeat steps 2 and 3 until two successive values for L agree. 

This model requires iteration for the basic variables B and L. In actual 
practice, many other parameters are to be specified. These can be com- 
puted directly once B and L have been determined . 

The solution here is interesting because the straightforward iteration 
procedure diverges. If equations (1) and (2) are stated as: 

B = f x (L) 

L = f 2 (B) 

then a successful iterative equation to replace equation (2) in step 3 above 
is: 

L i+1 = h J i 2{ B i + l) -h \ (3) 



where 

B 1+1 = f^) 

Using this iterative procedure, convergence is relatively slow. An 
attempt to improve the convergence using the Newton-Raphson technique 
gives a divergent iterative scheme. However, a modification of the pro- 
cedure shown does improve convergence. 

The improvement was made simply by changing the equation to read: 

L i+1 = L i - < f 2 ( B i + l> " L i> < 4 > 

This is treading on obviously dangerous ground, since arbitrary changes 
in the iterative equation are not easily justified. For this problem, how- 
ever, it does work. This is mentioned here to show that experimentation 
with techniques is necessary in some instances when other approaches 
fail. But keep in mind that it is important to have some way of judging 
the correctness of the results. 

Consideration of the following graph may assist in understanding the 
iterative process. Let 

x n+l = f (x n ) 
be the direct iteration equation, and represent all other schemes as 
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x n+1 = W(f (x n ), x n ) 

where the W represents the weighting function used in making the next 
estimate. 

In either case the evaluation of the function giving the next estimate might 
be plotted as in the graph. 

For a solution, the next estimate of x and the input value of x must be 
equal ( f (x) and x as shown in the graph). Assume an error, € , is made 
in the first estimate, x . The graph indicates the next estimate to be in 
error by an amount 6 2 where 6 2 > ^ 1 (* ne technique is diverging and 
6 3 is even greater than 6 2) . 

Given the initial estimate x n and the graph as shown, an obvious proce- 
dure for preventing the divergence is to make use of the difference 
( € 9 - 6 , ) in writing a new iteration equation — that is: 

x n + l =x n-< e 2- € l> 
and since ( 6 2 " ^ 1) = f ( x n^ " x n 
then 



x , _ = x - (f (x ) - x ) 
n+1 n v x n' n' 




ei (ea-eo (e 3 -e 2 ) 



This is recognized as equation (4) above for the cam problem. If 
( 6 2 - ^ 1) > ^ 1 this procedure overshoots the correct result and 
divergence may again result. Equation (3) in the cam problem used a 
weighting factor of 1/2 to prevent the overshoot. Choice of the "best" 
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weighting factor would require complete knowledge of a similar graph — 
which of course is not available without a solution. 

The intent here is to indicate the reasoning which may go into the selec- 
tion of iteration equations . 

The direct and Newton-Raphson iterative techniques are not the only pro- 
cedures available. For example, wide use has been made of a procedure 
credited to Mr. J. H. Wegstein of the National Bureau of Standards. The 
Wegstein method does not require evaluation of the first derivative and 
therefore is a powerful tool in cases where the first derivative is im- 
possible or difficult to evaluate. 

The rocker arm cam problem illustrates the use of experimentation to 
find a successful iterative scheme. No mention has been made of the use 
of an iterative technique to handle simultaneous algebraic equations. 

The problems discussed in this manual have been treated in the form 

x = f (x, y, z) 
Physical problems may frequently give rise to forms similar to: 

x = f-L (x, y, z) 

y = f 2 (x, y, z) 

As a matter of fact the rocker arm cam problem can be treated as a prob- 
lem giving rise to simultaneous algebraic equations . 

The Newton-Raphson and Wegstein methods may both be extended to han- 
dle simultaneous equations. 



Spring Problem 



The engineering design of springs makes an excellent computer applica- 
tion, especially from the economic viewpoint. It takes a great many 
springs to keep this mechanical world operating, and a less expensive or 
longer-lasting spring design can mean large returns dollar-wise. This 
example is concerned with the design of a particular type of spring, but 
the general method of computer solution would apply to almost any spring 
problem . 
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i.n D-L 



TD DT 



DD 



DD 



Cross Section of a 
Double Helical 
Spring 



A double helical spring is to be designed to support a given torque winding, 
such that the mean diameter and length of the spring will be within set 
limits — the shorter the length the better. The variables specified (sub- 
script 1 refers to inner spring, subscript 2 to outer spring) are: 

M-^ , M 2 = Winding torque 

S- , S„ = Maximum selected stress, or ultimate stress, of material 

T^ , T 2 = Maximum torsion angle 

K-. , K 2 = Required spacing between turns 

C = Required spacing between coils 

E = Young's Modulus 

The following are to be computed: 

h-i , h 9 = Dimensions shown in the diagram 



b 1 



Dimensions shown in the diagram 



D^ , D2 = Diameter of springs 

L^ , L,£ = Length of springs 

N-j_ , N£ = Number of turns in the springs 

Known relationships are: 
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(1) M = E " n T Gives winding torque for a given torsion 

6.6DN angle . 

Sbh2 

(2) M = — — Limits winding torque in terms of ultimate 

stress of the material. 

(3) L 2 = (N 2 + l)(b 2 + K 2 ) Outer length by geometry. 

(4) Di = D2-hj_-h2-2C Inner dimension by geometry. 

Note that for the outer spring five variables (h2, b2> D2, L2, N 2 ) must 
be determined but there are only three equations which apply. One 
approach to obtaining a solution is: 

(a) Set Do (because limit exists) 

N 2 (because this must be either an integer or 

an integer plus . 5 — successive guesses 
will be simpler) 

(b) Equate (1) and (2) and solve for h. 

= 1. 1D 2 S 2 N 2 
ET 2 

(c) Then from (2) obtain 

6M 9 
b 2 = H (6) 

s 2 h 2 2 

(d) and L 2 = (N 2 + l)(b 2 + K 2 ) (7) 

This provides the parameters for the outer spring. At this point it is 
necessary to make an additional guess to handle the inner spring. 

(e) SetN 1 

(f) Then solve 

D 1 = D 2 -h 1 -h 2 -2C 
and h x = l.lDjSxN! 



ET X 




simultaneously for h-, obtaining 




l.l Sl N l( D 2 -h 2 -2C) 




ETj + l.lS^ 


(8) 


(g) From (2): 




6M-L 

b 1 2 
1 S l*l 2 


(9) 
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(h) and finally: 

L 1 = (N 2 + l)(b 2 + K 2 ) 

Now testing is performed to see whether the springs are suitable; L 2 
must be equal to L, . 

According to the design criterion — that L 2 approximately equal L, — a 
decision can be made to accept or reject this design. In case of rejection 
the selection for D 2 , N 2 and N]_ can be modified and the indicated cal- 
culations repeated. 

A logic diagram of the procedure is shown in Figure 7. 





Set 
D 2 , N 2 , Nj 














+ i 






For Outer 
Spring Calculate 

h 2' b 2> L 2 










Select 

New 

D 2 , N 2 , N x 


1 


r 




For Inner 






Spring C 
hj, b, L 


alculate 
I 






A 




Satisfactory 
Design 



Figure 7. 



The question of how to modify the selection of D 2 , N 2 and N* must be con- 
sidered. The mathematical nature of the problem (nonlinear equations 
with an infinite number of solutions) does not permit use of the iterative 
techniques mentioned thus far. Two possible procedures for varying the 
selection may be called scattering and search. Both of these procedures 
assume the existence of a maximum or minimum for the function G which 
is used to denote "goodness" or "acceptability." In the spring problem 
the function is 

G= I L 2 - LJ =G (D 2 , N 2 , N x ) 
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Considering G as a function of only two variables 

G = G (D 2> N 2 ) 

allows visualization of the graph of G shown here. A corresponding three- 
dimensional picture for 

G (D 2 , N 2 , N x ) 

could also be visualized. 

The continuous curves represent contours on which G is everywhere equal. 
There may exist one or more minimums as shown. 

Scattering 

A set of selected values for D 2 and N 2 (in three dimensions, D 2 , N 2 , and 
N-, ) which create a grid covering the possible range in both D 2 and N 2 is 



T 



Allowable 

Range for 

D_ 



Allowable Range For 
Figure 8. Scattering 



O 1st trial points 
A 2nd trial points 




used to compute corresponding values for G. The grid is shown in Figure 
8 as a set ofO's. A smaller area is selected around the minimum G and 
the process is repeated, as indicated by the A's in the figure. 

This corresponds to the familiar trial-and-error approach except that 
many tries may be run at one time on a computer. Each solution may be 
made available allowing consideration of secondary design criteria. In 
the spring problem, a minimum | L 2 - Lu | is not sufficient for a good 
design; the ratios^! » h^and h_i are also important. It is quite probable 
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that the finally selected spring would have a small but not the smallest 



| L 2~ L i 
Search 



and other desirable features. 



One starting point in n dimensions is selected (shown in Figure 9 for the 
two-dimensional case). 




Figure 9. 
Search 



Start 



For each variable in turn: 

(a) The variable is changed by some small amount and G is recomputed. 

(b) If G decreases, (a) is repeated until such time as G begins to increase 
and then the variable is set to correspond to the smallest G. 

(c) If G increases for (a), an equal change is made to move in the opposite 
direction in a manner similar to (b) . 

Steps (a), (b) and (c) are repeated for each variable in turn, until such 
time as new changes within the practical bounds fail to decrease G sig- 
nificantly. 

This method moves more quickly to the solution and represents a more 
fully automated procedure . However , it demands that the behavior of the 
function (in this case G) be well understood. Discontinuities or lesser 
"minimums" can destroy the effectiveness of this method. 

The scattering approach can often be used to define the behavior of a 
criterion function such as G (D£, N£, N^) above, to the extent that the 
search technique may be used in succeeding computations . 

It is possible to add other selection criteria in the search technique. For 
example, the material cost in designing the above spring might be used 
such that the selected spring has the lowest material cost for those designs 
with a value of G within a specified range. 
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Pitman Arm Stress Analysis 



A problem often encountered in mechanical engineering is that of design- 
ing a member to withstand given tensile and shear stresses. This partic- 
ular example concerns a pitman arm (see Figure 10) as a part of an auto 
steering mechanism. The drawing shows a pitman arm with the applied 




SECTION AA 



Figure 10. Pitman Arm 

force F which gives rise to the various stresses in the arm. 

The usual approach to this problem is for the engineer to specify dimen- 
sions which seem appropriate and then to check the design with calcula- 
tions. This is done by choosing cross sections of the design at intervals 
along the arm — computing the stresses at each cross section and check- 
ing that these stresses are less than the failure stress for the material 
used. 

For Section AA shown, the computations are: 

(a) For bending stress 
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s b = £k (i) 

Z x 

For torsional shear stress 

h = ^ < 2 > 

z 
p 

where 2? x and Z are called section moduli and depend upon the moments 
of inertia for the member. For the elliptical case the moduli are 

Z x = f 2 bh 2 (3) 

z p ■ fe bfa2 < 4 > 

(b) The total normal tensile stress and the total shear stress can be com- 
puted in terms of the bending stress and the torsional shear stress. 

For normal tensile stress 



A. T 2 2 1/21 

5n = 2 l-Sb +(S b + 4S t ) J 



(5) 



For shear stress 



12 2 1/2 

S s=^ S b +4S t> < 6 > 

Once these total or maximum stresses have been calculated, they are 
compared with acceptable values for the material to check whether a 
feasible design has been arrived at. 

An approach which may be of more use to the engineer — and which is 
facilitated by the use of computers — is to start with the maximum allow- 
able stresses and from this to obtain the dimensions required for the 
member to support these stresses. 

For example, assume that the pitman arm shown is to be elliptical in 
cross section, with h=2b, and that values have been assigned to S n and 
S g . The problem now consists of determining values for b and h in terms 
of these values. 

Equations (5) and (6) can be manipulated to give S b and S t in terms of S n 
and S s : 

Sb = 2(S n - S s ) (7) 

r 2 2l i/2 

S t = L S s - <S n -S s > J (8) 

Equations (1) to (4) may be used finally to restrict b and h in terms of 
S b and S^: 

w uu 2 F1 
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and 



" 9 



^bh-=Fli 



and 



b 6 -_8 Fl 2 = 

IT S^" 

b 3 - 4 Fli = 



or, since h = 2b, 

(9) 
(10) 



7T S 



This gives two cubics which restrict b. The proper design must have a 
value of b which is at least as large as the larger of the two roots obtained 
from equations (9) and (10). 

Now the proper dimensions can be computed for the arm as we move along 
the arm (that is, change 1^ and 1 2 ). Furthermore, if the curve used for 
the arm can be expressed mathematically, then 1^ and 1 2 can be auto- 
matically computed as the arm is traversed toward or away from the axis 
of rotation. If this is impossible, then 1, and 1„ must be measured for 
various cross sections of the arm and the corresponding dimensions com- 
puted for the cross sections. 

The point of this example is that very often a computer approach should 
start with a reanalysis of what is to be done, considering the possibilities 
of improving the usefulness of the results. 



Automobile Suspension 



The following three-dimensional problem occurs in the engineering of auto 
suspension systems; arms of length d]_, d 2 and d 3 are fixed at points in 
space A,B and C respectively. 



(x,y,z) 
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The point of juncture of the three rigid arms is found to be the intersection 
of three spheres with centers at A,B and C and radii d^, d 2 and d3 respec- 
tively. 

Expressing this in equation form gives the following three equations: 

(x - x x ) + (y - y x ) + (z - z x ) = d x (1) 

(x - x 2 ) 2 + (y - y 2 ) 2 + (z - z 2 ) 2 = d 2 2 (2) 

(x - x 3 ) 2 + (y - y 3 ) 2 + (z - z 3 ) 2 = d 3 2 (3) 

This set of equations appears difficult to solve for x, y and z, and iteration 
has been used in many instances to find a solution. Iteration in this prob- 
lem is time-consuming and unreliable (convergence is not always possible) . 

If the following algebraic steps are performed, an explicit relationship for 
x, y and z is found: 

(a) Perform indicated squaring in equations (1), (2) and (3). 

(b) Subtract equation (2) from (1) and similarly equation (3) from (2) . 
This will eliminate terms in x 2 , y 2 , and z 2 giving two linear equations 
in x, y and z. 

(c) Solve these two linear equations algebraically to obtain x and y as 
linear functions of z. 

(d) Substitute these functions for x and y into equation (1) and solve the 
resulting quadratic in z. 

This will then give equations of the form: 

9 -i 

a = -b + (b^ - 4ac) 2 

2a 
x = e + f z 

y = g + hz 

where a,b,c,e,f,g and h are intermediate substitutions for terms made 
up of the known quantities xj, yj, Zj and dj. 

This gives two complete sets of answers: the correct physical answer and 
its mirror image as illustrated by the dotted lines. The selection of the 
proper set of results can be accomplished by programming the logical 
considerations of the physical system . 

A computer is convenient for the evaluation of the lengthy equations in- 
volved. If the problem is to accomplish this evaluation for many points 
in the motion of one of the points (say point C, as indicated by the arrow), 
the use of a computer becomes mandatory. The orientation of a suspen- 
sion system must be known for the many possibilities in the auto's motion. 
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This problem illustrates an important class of mechanical engineering 
problem — that of determining the orientation of physical bodies by speci- 
fying the coordinates of points in the body. 



Spring Stress Analysis Problem 



The following illustrates the use of the same algebraic approach in another 
problem which also involves the use of the summation of forces so familiar 
in stress analysis. Find the necessary tension, S, in the spring for the 
support of W at various values of the angle -6- . 




The analysis might proceed as follows: 
To find point A 



x l = ^1 cos "®" 



y 1 = d 1 sin 



To find point C 



. ,2 . ,2,2 
(x-x-l) + (y-y-L) = d 3 

(x-x 2 ) 2 +(y-y 2 ) 2 =d 2 2 



(1) 
(2) 

(3) 
(4) 



Following the procedure as outlined in the previous problem, expand (3) 
and (4) to give: 



x - 2xx-L +x 1 + y - 2yy x + y 1 = d 3 
x 2 - 2xx 2 + x 2 2 + y 2 - 2yy 2 + y 2 2 = d 2 2 
Subtract (6) from (5) to give: 



(5) 
(6) 
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2x (x 2 -x 1 ) + (x ± 2 - x 2 2 ) + 2y (y 2 -yi) + (yi 2 -y 2 2 ) = (d 3 2 -d 2 2 ) 
or x = Ey+F (7) 



Where 



E = 



fri-y2> 
(x 2 - Xl ) 



and F = [(d 3 2 -d 2 2 ) - (y 1 2 -y 2 2 ) - (x^x/)] /2 (x 2 - Xl ) 

Substituting from (7) into (5) gives 

E 2 y 2 + 2EFy + F 2 -2x 1 Ey-2x 1 F + x 1 2 +y 2 -2yy 1 +y ][ 2 = dg 2 or 
(1+E 2 )y 2 + (2EF-2x 1 E-2y 1 )y + (F^x^+x^+y^-dg 2 ) = 



or finally y = - Q + Vq 2 -4PR 



2P 



where 



P = (1+E 2 ) 

Q = (2EF - 2x x E-2y 1 ) 

R = (F 2 -2x 1 F+x 1 2 -y 1 2 -d 3 2 ) 
and x = Ey+F 
To find angle (J) 

x 3 = ^1 + d 4^ cos "®" 
y 3 = (dj + d 4 ) sin -9- 

and(() =tan" (y-y 3 ) 

The determination of (|) allows the use of the force diagram and the equi- 
librium principle, which states that the sum of forces acting in both the x 
and y directions must be equal to zero, to determine the necessary tension, 
S, in the spring. 




*W 
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V F = T cos-9- - S sin ( (J) - 90°) = 
JF =T sin -6- + S cos ( (J) - 90°) - W = 



or 



T = 



S sin ( $ - 90°) 



cos-O- 

S = W/ [cos ( (j) -90°) + tan -9- sin ( (J) -90°)] 

The usefulness in the design of a program to evaluate the tension in the 
spring for any configuration, motion of-G- , or load W is obvious. In this 
case the evaluation depends upon locating points in a plane. Traditionally 
this type of determination has been made on a drafting board. As in this 
case, where motion or many possible designs are to be evaluated, cal- 
culations on a computer can greatly reduce the time invested. The improve- 
ment in accuracy can be even more important. For example, a one-degree 
error in the steering geometry of an auto can completely change its steering 
characteristics . 



Vehicle Simulation Model 



It would be unfair to close this brief discussion of a few engineering 
mathematical models without mention of at least one of the larger models 
which have been constructed. 

Several auto manufacturers have succeeded in representing the operation 
of cars and trucks by a computer model. By operating the model in the 
computer, such items of performance as acceleration, fuel economy, and 
riding characteristics can be determined for a proposed model before the 
model is built . 

The procedure for accomplishing this introduces a very interesting con- 
cept. Normally the solution to a set of equations is obtained by solving 
for one point in time — a static situation. The performance of a car, 
however, is a dynamic situation. Beyond this, equations to describe all 
the myriad operations in moving a car simply do not exist. Rather, there 
exist combinations of theoretical equations, empirical equations and tab- 
ular data to describe the functioning of various parts of a car. The heart 
of a car simulation model is the logic which ties these pieces of theory and 
data together under one unifying variable. Different approaches choose 
different unifying variables; for illustration, however, suppose the variable 
chosen is time. 

Typically, the motion of the car might be studied time-wise for a set level 
of fuel input. For the first unit of time the information representing the 
relationships between parts of the car is interrogated for the motion to be 
expected, working from the engine back through the transmission and 
finally to the wheels. Once the motion for the first unit of time is deter- 
mined, the time is incremented and the operation repeated. 
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Another approach might be essentially the reverse of this: Starting with 
the desired velocity of the vehicle, work back through the wheels, dif- 
ferential and transmission to determine what the engine speed must be to 
give this velocity. However, time is again the unifying factor in that con- 
ditions for the previous time interval help determine the required motor 
speed for a particular velocity. 

A look at the basic relationships in this example may be helpful. 

Assume we want to know the linear acceleration of the vehicle at full 
throttle as a function of time . 




Starting with a given engine speed N^ and linear velocity V, the first set 
of relationships is concerned with the engine, A. 

(a) The first step is to use a table to determine the output torque of the 
engine . 

T i = f i< N i> 

The torque delivered to the transmission may likewise be a tabular 
function of N^: 

T 3 = f 2 ( Nl ) 

as is the torque consumed in friction of the engine parts: 

T 2 = f 3 < N 1> 

At this point the acceleration of the engine for this time interval may 
be determined as 



©< 



T r T 2- T 3 



m j 

where I = total engine inertia 

(b) Drive shaft speed, N 2 > is computed by 
N 2 = V_ G 



where G = axle ratio 
and L = rolling radius. 
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The transmission ratio, TR, is a function of N„/N.. 

TR = f (Ng/N^ 

This relationship is best represented by a table of TR vs Ng/Nj which 
has been determined experimentally for the particular transmission 
to be used. 

(c) Moving on to the drive shaft, the torque delivered to the drive shaft 
may be computed as 

T 4 = T 3 . TR 

(d) Finally, arriving at the rear wheels, the force exerted on the wheels 
due to the delivered torque may be given as 

F =T < T 4 - T 5> 

where T 5 = torque loss due to friction and acceleration of drive line 
parts . 

Then the acceleration of the car, a, is given by Newton's second law 
of motion 

a = F -R 

M 

where M = mass of the car and 
R = rolling resistance. 

The rolling resistance is generally given as an empirical function 
such as 

R = aM + b AV 2 

where A = frontal area of vehicle, and a and b are experimentally 
determined constants (see section on empirical relationships) . 

At this point the time is incremented and a new engine speed N 1 is com- 
puted as 

*-m = ^m + A * m+1 



^i> 4. i = (N iL + ^\ °< At 

1 m + 1 -Lm \ m 



m 



The above procedure A,B, C and D may be repeated for as many time 
increments as desired. 

This logical procedure may be represented by the flow chart in Figure 11, 
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Figure 11. 



Figure 12. 



The flow chart in Figure 12 is included to illustrate how essentially the 
same mathematical model may be used in reverse to compute fuel economy. 
In this case the desired motion as a function of time is specified. Cal- 
culations determine the engine speed required to give the desired motion 
for each time interval. Finally the fuel consumed at a given engine speed 
is known and therefore the sum of the fuel to be consumed for all time 
intervals may be calculated. 

Many simplifications of vehicle simulation have been made in this descrip- 
tion. 

The major objective in presenting a description of vehicle performance in 
this chapter is to illustrate the use of the three types of information 
(theoretical equations, empirical equations, and tabular test data) under 
a logic format to obtain a computer simulation of the physical object. 

The next chapter describes the use of the computer as a tool in finding 
empirical relationships . Actual procedures for handling of tabular data 
in a computer program are also discussed. 
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Exercises 



1. Perform the algebra required to obtain equation 9 in the solenoid de- 
sign problem , using equations 1 to 8 . 

2. Perform the algebra indicated in the suspension analysis problem and 
either show the complete equations for x, y and z or show the equa- 
tions which define the intermediate variables a, b, c, r, s, t and u. 

*3. Numerical techniques are often used to solve simultaneous differential 
equations with nonlinear coefficients . Prepare a logic flow chart to be 
used to program the solution of: 

d2 y , dx 

J +o( — +xy = A 



dt2 dt 

dx+dy +ot =B 
dt dt 

where A and B are known constants and oC = t + sin oC t, given t QJ x , 
y and fdy\ . 
W 

Hint: Use Euler's approximations as described in problems I and J 
of Section IV for the differential equations and Newton's procedure 
for evaluating &C at each interval in t. 
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SECTION VI: SELECTED MATHEMATICAL TOPICS 



The problems discussed in this section have been chosen to answer some 
of the legitimate questions that many practicing engineers might well ask 
about terms often heard in computer groups. Such things as matrix inver- 
sion, relaxation techniques and eigenvalues are out of the realm of expe- 
rience of many engineers. It is felt that some simple examples illustrat- 
ing these terms may lead to eventual recognition of real problems which 
can be tackled using the techniques described. 

Matrices 

A matrix is an array of numbers . Consideration here will be restricted 
to (1) the rectangular matrix, which may have the form: 



A = 



11 d 12 



a 21 a 22 



a ml °m2 



L ln 



2n 



mn 



where m is the number of rows and n the number of columns; and (2) to the 
column (column vector) matrix, with the form: 



B = 



m 



Matrix Addition and Multiplication 



One of the most common applications for matrix notation is in trans- 
formations. A simple transformation can be expressed as the addition 
of two single-column matrices. The figure below illustrates transformation 
(translation) of the coordinates for point P from the x, y, z coordinates to 
thex', y', z' coordinates. 



P(x, y, 2) 
(x', y', z-) 




82 



If the coordinates of P in unprimed coordinate systems are represented by 
the column matrix. . . . 



A = 



. . . the coordinates of O'are represented by . . 



B = 



x o 

y 

z o 



. and the coordinates of P in primed systems by .... 



C = 



x' 

y' 

z 1 



. . . transformation (shift in space) from unprimed to primed systems can 
be expressed as .... 

C = A - B 



or 

The inverse transformation is .... 

A = C + B 
x 1 +x. 



- 




r~ ~] 


x' 




x-x Q 


y' 


= 


y-y 


z' 




_z - z _ 



o 

y' + y 

z' +z 



The above example illustrates a matrix addition. 

Another type of common point translation is rotation about an axis . 
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The illustrated rotation of the point P about the Z axis through the angle -0- 
can be written .... 

x' = cos x + sin y 
y 1 = -sin0 x + cos y 
z 1 = z 



. . . or in a 


matrix notation 


A = BC 


where . . . 


A = 


x' " 

y 

z 1 


B = 



cos 


sin 





-sin 


cos 











1 



c = 



which illustrates the use of matrix multiplication. 

Frequently it becomes necessary to transform from an unprimed coordinate 
system to a primed coordinate system which is both linearly and rotationally 
different. This requires both matrix addition and multiplication. 



• P(x, y, z) 

(x 1 , y, z') 




Given the coordinates (x, y, z) of P, the transformation to primed coordi- 
nates is .... 



x 1 = X i x + ui y +v^z + x Q 
y 1 = \ 2 x + u 2 y +v 2 z + y Q 
z' = X 3 x + u 3 y + iTg z + z Q 
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Transpose 



where (Xi ,ui,vj) , (X2»u 2 ,^ anc * ^3 ' u 3 ,I/ 3) are the direction cosines 
of the x\y\ and z' axes respectively, in the X, Y, Z coordinate system. 
In matrix notation this is 



A = BC 



D 



The above examples have been included because transformations are very 
useful in handling many problems of motion of physical systems. The most 
common form of mechanical motion is rotation about an axis. 

Computer programming of a complex geometrical system is greatly simpli- 
fied if the geometry and motion of the system are expressed in the most 
general mathematical terms. The use of trigonometry, because it forces 
a "one-position-only" or static description of such a system, may seriously 
retard a satisfactory computer handling of the problem. 



The transpose of a matrix is obtained by interchanging rows and columns 
in the matrix. 

The transpose of . . . 



B = 



Xi 


ui 


vi 


^2 


u 2 


v 2 


*3 


u 3 


v 3 



from the previous topic is . . . 



B' = 



u l 



^2 



u 3 
v 3 



Then the inverse transformation from primed to unprimed coordinate 
system for the previous figure can be expressed as . . . 

C = B' (A -D) 

where A, B, C and D have the same meaning as defined above. 

Matrix Inversion 

A major application of computers is in handling the solution of large sets 
of simultaneous equations which may occur in such engineering areas as 
stress analysis, statistical least squares and circuit analysis. 



E, -=k 



•kz.iT R 3 1 yr ^z. E ' 
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For the simple circuit shown, applying Kirchhoff's law gives a set of linear 
equations in i and i : 



(R-j: + R 3 ) i-^ - R 3 i 2 - Ej 
-Rgi! + (R 2 +R 3 ) i 2 = -E 2 



This set of two equations in two unknowns does not require sophisticated 
handling. However, consider a slightly more complex circuit in which 
the values of the resistances are known and the currents are to be deter- 
mined. Again, Kirchhoff's law may be used to establish a set of linear 
equations: 

I x + I 2 + I 3 = Ai 



*8 + *9 + ^0 



- I 1 +I 4" I 6 



-13 + I5 " *9 



W 1 !*) 



= A, 



= 



1 amp 



= 



-R7I7 + RqIq " R 10 I 10 ~ ° 



-R 5 I 5 + R 8 I 8 - R 9 I 9 - 




R 2 X 2 " R 3 I 3 " R 5 I 5 



-R-iI-| + R9I9 ~ R 4^4 ~ ^ 



" R 4*4 ~ R 6^6 + R 7*7 ~ ® 



Since the solution of this set of linear nonhomogeneous equations requires 
a considerable expense of time and effort, the problem is well suited for 
a computer. 

Matrix notation is particularly useful when dealing with large sets of linear 
equations. For example, the above set of equations may be represented 
by 

AX= B 

The solution of such a set of equations may be obtained by one of several 
exact methods or by one of many iterative methods . 
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Although both types are used in computer programming, the exact is pre- 
ferred for its accuracy and reliability. However, the size of the equation 
set may prevent the use of an exact method. The iterative method requires 
considerably less computer memory for the handling of simultaneous 
equations. 

In matrix notation, the solution of the equation 
AX= B 

is equivalent to finding the inverse A of the matrix A. One of the 
properties of A is 

A"* 1 A = I 



where I is the identity matrix with the form: 



10... 

01... 

. 1. . 

1. 

1 



where the diagonal elements have a value of 1 and all other elements are 
zero. 

The I matrix has as one of its properties that premultiplication by the I 
matrix leaves any matrix undisturbed as 



IX = X 



Therefore, if both sides of 





AX= B 




are multiplied by 


A" 1 






A _1 AX = A" 


-1b 


the result is 


X = A _1 B 
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Thus , one method which may be used in finding the solution of simultaneous 
equations is calculation of the inverse of the coefficient matrix. 

An exact method for matrix inversion is shown below. The problem on 
eigenvalues discussed later in this section demonstrates an iterative 
procedure for the solution of simultaneous equations. 



Inversion by Elimination 

Given the coefficient matrix: 



a ll a 12 a 13 • • • a ln 
a 21 a 22 a 23 • • • a 2n 



a nl a n2 a n3 ' ' * a nn ° 



□□□...□ 



add a column — a unit vector — which contains a 1 in the first row and 
zeros elsewhere. At the same time, add a row — called the pivot row - 
denoted by CD 's. Then perform the following computations to arrive at 
a new array: 

1. For the pivot row elements: 



P.J 



l l, 3+1 
a l.l 



j = 1,2, ...,n 



2. For all other elements, compute a new value: 



a' =a -(a ) (a ) i=l,2 n 

i,j i,j+l i,l p,j j=1>2 ,...,n 



3. As a result of step 2, all the new elements of row 1 are zero. This 
row is dropped, and the remaining n rows renumbered 1 through n. 
Thus, for the last row: 



a' = a 
n,j p,j 



j = 1,2 n 



4. Add a new unit vector and pivot row and repeat steps 1, 2 and 3 a total 
of n times. The resulting array is the inverse of the original matrix. 

Given a set of simultaneous equations: 
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a ll x l + a 12 x 2 + 



+ a ln x n = b x 



a nl x l + a n2*2 + • • • + a nn x n = b r 



then the solution can be obtained directly by starting with an n+1 by n array 
in which the original matrix is augmented by the b vector. If values of x 
are required for more than one set of b values , the additional b vectors 
can be incorporated in the original array as shown: 



a-M a-io . . . a-i n b-M b 



L ll d 12 



In "11 u 12 



a 21 a 22 ' * * a 2n b 21 b 22 



a nl a n 2 • • • ann b n i b n 2 



Subjecting this n by n + 2 matrix to the four- step procedure above a total 
of n times gives the array: 



(- 



x n x 12 a^ aj_ 2 . . . a.\ 



x_ x 



21 22 21 22 * " * 2n 



xx a' a' 
v nl n2 nl n2 




where the a! . are the elements of the inverse of the original coefficient 

matrix, and the xj s j are the solutions for each of the two b vectors of the 
matrix equation 



AX =B 



Eigenvalue Problem 



This problem is discussed to show, with a simple example, an application 
of computers in handling sets of differential equations . It illustrates the 
meaning of eigenvalues for a set of differential equations which describe 
the motion of a mechanical system. 



A 



*i 



/ y — jimsub—C^) — <ms&r-Q\-<Mmsu--CX-4s 



x 2 



WUUULr— 



\ 

', 



m^m^—m- 



k 1 =k 2 =k 3 =k 4 
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Given the spring and mass system shown, assume that the normal modes 
of oscillation of the system are to be determined. 

The differential equations of motion can be written: 

m-jx^ = - kxj + k(x£ - X]j = k(x£ - 2x-^) 

m 2 x 2 = " k ( x 2 " x l) ~ k ( x 2 ~ x 3) = " k ( 2x 2 " x l " x 3) 

m 3 x 3 = ~ ^3 + k ( x 2 ~ x 3^ = k ( x 2 ~ 2x 3^ ^ 

A procedure for solving these differential equations is to assume solutions 
of the form 



x = Ae iaft , x 9 = Be ia/t , Xo = Ce iujt 



(2) 



Substituting these into equations (1) , performing the indicated differentia- 
tions and rearranging terms gives: 



(3) 



(4 


- u?) A - 


h b 


= 




-1 A + <4 


- 2 )B-|C 


= 




- 


k- r» + /2k _ 
m B + (m" 


(jj 2 )C = 



which can be written in matrix form as 
/7 2 m-m <A A 2 



=K + 2k _ k_ 
m m m 



0-JL + 2k 
m m 




or more simply, 

(D- "M)X = 



= 



(5) 



(4) 



For A, B, C to satisfy equations (3), the determinant of the coefficient 
matrix must vanish. This is equivalent to the statement 



(D- M)=0 



(6) 



Evaluating the determinant for equations (3) and equating it to zero gives 
a polynomial — called the characteristic equation — in ou 2 : 



x m 
with roots 



&■ - oa 2 ) 3 - 4 (4 " «> 2 ) = 



m x m 



(7) 



2_2k 

1 m 
2_ 2k 

2 ~ m 



w 



2 2k /— k 
a) . = — + V2 — 
m - 1 



m 



(8) 
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The values of oj which satisfy these equations are called eigenvalues. In 
general, values of >\ which" satisfy (5) are eigenvalues. The vector X is 
called the eigenvector. 

It can be seen that for larger sets of differential equations the solution of 
the characteristic equation requires some means for computing the roots 
of a polynomial. 

For this problem the eigenvalues give the natural modes of vibration for 
the mass and spring system. This calculation of eigenvalues for differen- 
tial equations is termed frequency analysis . In an earlier section attention 
was given to the solution of differential equations. In motion problems 
this amounts to determining the displacements as a function of time and 
is termed an amplitude analysis . Both frequency analysis and amplitude 
analysis are important computer applications. 

To illustrate a computer solution to a frequency analysis , consider the 
spring and mass system shown in the diagram. 




The differential equations which describe this system are as follows: 
d x- 



dt 



2T + a ll X l + a 12*2 + a 13^1 + *14% = ° 



d xo 

+a x + a x +a -9- + a_ -9- =0 

,.2 21 1 22 2 23 1 24 2 

dt 



2 

— - + ag^ + a 32 x 2 + a 33 ± + a 34 2 
dt 2 



d 2 -e , 

dt 2 



2 + a 41 X l + a 42 X 2 + a 43 e i +a 44^2 = 



where the a. . are dependent on the spring constants and the masses, 
Assume solutions of the form: 



x, = X_ COS (tit 

1 10 



9- = 9- cos CD t 
1 10 
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x 2 =x 2() cos wt 



&2 = ^20 C0S <* 



where x 1Q , x 2 q, -9- 10 and -©-20 are ^ e initial displacements. The appro- 
priate differentiations and substitutions give a homogeneous set of linear 
equations of the form: 



a ll a 12 a 13 a 14 



a 21 a 22 a 23 a 24 



a 31 a 32 a 33 a 34 



a 41 a 42 a 43 a 44 





X 10 




x 20 




1O 




" 9 '20 



OJ 



10 



k 20 



10 



r 20 



The iterative procedure for the solution of these equations involves the 
following steps: 

1. Initially, let x 10 = x 2Q = -9- 10 = -9- 20 = 1.0 and evaluate the left-hand 
side for new values of u^x^q, aJ 2 x 2 o> tL)2 '®'io anc * Oj2 "®'20' 

That is: 

2 
oix = a x +ax + a -a- + a -a 
10 11 10 12 20 13 10 14 20 

and so on. 



2. "Normalize" for new guesses at x 1 , x ? , -a- , and -B-^n k y setting 



X 10 = 1 



cu 2 -e- 



1P_ 



<A 



10 



u) u x 



X 20 = 



20 



10 



2 

Gl> X 



10 



2„ 
uj -9- 



20 



k 20 2 
(JJ x 



10 



3. Repeat steps 1 and 2 until such time as successive values of x^q, x _, 
•0- , and -e- are very close. At this time, convergence has occurred 

and o> 2 can be computed. 

2 2 

4. Upon convergence uj can be evaluated from the given value of oi x-^; 



say 



U x 10 = K 
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The following is the iterative solution for a> 2 of a sample set of coefficients: 
Particular equation set 



+929 


-332 


+3,590 


-748 


+748 


-8,090 


+4.17 


-4.17 


+35,400 


-14.31 


+14.31 


-40,300 



-3590 




x l 




x l 


+8090 
-11,760 




x 2 
*1 


= u, 2 


x 2 
*1 


+40,300 




5 2 _ 




j>2_ 



Iterative solution 



Iteration 


x l 


x 2 


*1 


*2 


1 




1.0 


1.0 


1.0 


1.0 


2 




597.0 


0.0 


23640.0 


0.0 


3 




143085.0 


-321095.0 


1401773.0 


-1595813.0 


4 




76883.0 


-171908.0 


477975.0 


-844314.0 


5 




63414.0 


-141558.0 


349238.0 


-693154.0 


. 6 




60681.0 


-135398.0 


323511.0 


-662485.0 


7 




60002.0 


-133869.0 


317130.0 


-654870.0 


8 




59825.0 


-133468.0 


315460.0 


-652876.0 


9 


a, 2 


59777.0 
= 59777.0 


-133362.0 


315017.0 


-652347.0 



Knowing that x-^q is to be set to 1, then 

cu 2 = K 

To clarify what has been done here, remember that the set of equations 
has no constant term — it is a homogeneous set of equations. Essentially 
this means that there are an infinite number of solutions which satisfy the 
equations . This is certainly reasonable when the physical system under 
consideration is examined. In a vibration problem of this kind the initial 
displacements x-^q, X2q, 
ferent values . 



■0-20' an( ^ "®"20 mus t k e expected to take on dif- 



In the above procedure a value of x-. Q was selected which then fixed the 
values of the other variables X2Q, "^io> "®20 anc * a H° we dcu to be deter- 
mined. 

A very similar iterative scheme (without the normalization step) can be 
used to solve sets of nonhomogeneous linear equations. 



Partial Differential Equations 



Many engineering problems involve the handling of partial differential 
equations. Three classes have been distinguished as: 

1. Elliptical equations (describing potential fields) 

2 
V (|) = g (x, y, z) 
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2. Parabolic equations (describing heat flow and diffusion) 

V 2 (t) =kii 
it 

3. Hyperbolic equations (describing wave action) 

* 2 <D 



V 2 0=± 



it 2 



where V is the Laplacian operator in rectilinear coordinates: 
2 A a 2fi \2<t) ^20 
ix" i v z $z ? 



v 2 (J)= ^2O_ + _i20 + 
- 2 >Y 2 



A basic approach to handling partial differential equations when describing 
a particular material or space is to create a grid of points covering the 
space as shown for a two-dimensional space. 



Then, at any point the 
first derivative with respect 
to x can be approximated in 
two ways: 



m 



Ax 



or 



(1) 



H\ „<t>o _ 4>1 



m 



Ax 



Y 



42 
AX 




44 



The second derivative can be approximated as 



*x 2 



Ax 



(Ax)' 



(2) 



Derivatives in the y direction can be obtained in the same way. Using 
this procedure, many partial differential equations may be reduced to dif- 
ference equations amenable to computer handling. 



Potential Problem 



Apply this to the two-dimensional problem of finding the potential dis- 
tribution in a square whose sides are maintained at voltages 

(V i)b , (V 2 ) b , (V 3 ) b , and (V 4 ) b 

If there is no charge within the square , the potential distribution is de- 
fined by the Laplace equation: 



* 2 V + 3 2 V 



ix2 



*y 2 







(1) 
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< V 3>B 



(V«) 



l'B 



(V 2»B 



< V 4> B 



?! 



Setting up a square grid system to cover the square, for the general 
point A, the approximations for the partial derivatives become: 



<> v ) 

ix 


1 


4 o 
h 


> 2 v 

>x 2 


*~> 


V 4 + V 3 - 2V 
h 2 


Similarly, 


for the 


y dimension: 


i 2 v 


~ V 2 


+ V, - 2V 



U 2 



_ V o" V 3 



iy 



Then (1) becomes: 



V l + V 2 + V 3 +V 4- 4V o =0 



(2) 



This is the basic relaxation equation. It is applied in the following way: 

1. A first guess at the potential of each point on the grid is made on the 
basis of the known boundary conditions . 

2. Moving systematically through the points of the grid, the quantity 
called the residual is computed for each point and stored. The residual 
is given by: 
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IBM 



Fm *7*-n27-i 



FORTRAN CODING FORM 





Punching Instructions 


Poge J_ of ^y 


Progrom rWtf - D/ME/VS/^Jtl. fi£Ljn/?r/0Af 


Graphic 
















Cord Form # * 


Identification 
i .. i .... 1 


Programmer iDott 


Punch 


















T] 00 
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NUMli> 

1 5 


5 FORTRAN STATEMENT 
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PJMEfiSItti V(3 0< t 20) } A(20 ti 20y 




Pt.aD, 10 ,$,»/ , J>£,1 ; ( ( V,( L t J),.I=l . 


M).J* 


l.N) 














L=M-A 




















Ll*Hrl. ....... . 




















Z><* fl? 1-1,, M ........ 




















*.#. .f.fr j-=i. f v , 
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R(l r J)-0., 




















M 3 J*Z,J. ............. 




















Wf .3 /=^^i ............ 


















3 


A(7. f J).*.y.(j. t .J.-i t )+y.(z, t .J.+i) 


+V(l- 


X } J)f 


V(l+1 


,J)-4 


.P*V.(j,.J.). 








4 


K*2. ................ 




















/>* f, z=*,,z. , 




















0/ f. J=.2 V L1 




















If(/J,BSf(A(l jS),)- DfL )9,9, 


8 
















2 


K.»±. ........... . . 




















RP£L*k(x. f J.)Jit>..i0 ........ 




















V.(T t J.)=.*.U,J)*JlP£L, 




















P(l,J) = Q.,' ' . 




















ft/2 ',,./- J )^I^Jrl)+HDBL 










.... 1 ... . 










R (l r J-hl) **(I. ',.1+1.) +XJEL 










.• 










lilrttJfrtfXr^Jl+APZL , 


.... 


. . . . 






. . 


1 . . . . 


1 . . . . 





* A standard cud form, IBM electro 888157, is available for punching source statements from this form. 
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Punching Instructions 
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pogrom f w j _j, IM£ H S j $tf/)L RELAX flTI0M 
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Punch 
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FORTRAN STATEMENT 
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* (1*1 , J)=* (I+L,J) +*J>£.L. 



C4NTJNUE 



C,§ T4 (+ f ll) t K 



JUL 



P*IA/T ZO tt ({V(I, r J) 
STSP, 919 ■ 

f fVrl I . , ■ . 



=/,*; 



,j=*, 



Ml 



Figure 13. 
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R = V ± + V 2 + V 3 + V 4 - 4V Q 

Initially, equation (2) will not be satisfied since the potentials are 
guesses. 

3 . Again moving systematically and considering each point not on the 
boundary, the potential is adjusted to make the residual for the point 
equal to zero by applying the following equation: 

new V Q = V Q + R Q /4 

4. This affects the residuals of the surrounding points, so they are ad- 
justed by: 

new Rj = Rj + R Q /4 

5. Steps 3 and 4 are repeated until no residual is found whose absolute 
value is greater than some predetermined limit of accuracy. At this 
time the relaxation equation is satisfied and the potential distribution 
is known. 

It is possible to quickly write a FORTRAN program to do the necessary 
computation. For this problem let 

M = Number of points in grid on x axis (200 max.) 

N = Number of points in grid on y axis (200 max. ) 

V (I, J) = Potential at points on grid (initial guesses) plus boundary 
values 

R (I, J) = Associated residual 

DEL = Limit of accuracy desired. 

The resulting FORTRAN program is shown in Figure 13. 

With this example in mind, it is of value to consider how a more complex 
problem of this type might be handled. 



Temperature Distribution Problem 



For a cylinder made up of layers of material of different conductivities, 
find the steady-state temperature distribution. The outer surface tem- 
perature (boundary temperature) is known. Here the use of cylindrical 
coordinates facilitates analysis. In cylindrical coordinates the heat flow 
equation is: 

_1_ _J_ / JT_\ + _1_ _YT_ + a 2 T = k a T (1) 
r Jr I k I r 2 ^ ({) 2 ^ z 2 ^t 

NT 

For the steady state, where = 0, then 

it 
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1L 
3r 




^ 2 T 



*f 



S 2 T 
^2 



= 



(2) 



Since this deals with a nonhomogeneous cylinder, any approximation to 
the partial differential equations must show a dependence on the con- 
ductivity of the material. A suitable approximation may be derived by 
reference to the accompanying diagram . 



Heat 





T i + 2 






T b -X 


• 
k i+2 






t i+i L 


y^f^ 




T 
i 


• 


• 




• 


k i + i 


k 




k. 

1 




T i + 3 








• 








k i + 3 








« h 







First assume the boundary between materials to be midway between two 
points on the grid and introduce a boundary temperature 1\ . Then the 
quantity of heat which would pass from one square and (assuming incorrect 
initial assignment of temperatures) the quantity received by the second 
must be equal; or 

k i < T i " T b> = k < T b - T > 



or 



rp _ fc Ti + kT 

1 b ~ x 1 



k +kj 
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Then approximate 



or 



^r 1 



Yr\ 

^r 1 



= T b -T 
h/2 



1^ +k 



h/2 



An equivalent approximation redefining T^ may be made from the opposite 
side to give: 

/JT\ J*i + 1 \ /"T i+1 + T\ 



/ \ k i+l +k / 



h/2 



The final approximation to the first derivative may be taken as: 



It 

^r 






+f m* yW" 1 



-T i+1+ T^ 



h + l +k J 



ki 



T--T 
l 



The second partial may be taken as: 

} r 2 ^r ^ ^r ) 1 \^ +1+ k)[ h 2 /2 j \^)\^/ 2 
For uniform boundary temperature: 
^ 2 T 



}<D 2 







and *T ^ / V3_ \ / T. +3 -T \ + / k i+2 

J * WW v h 2 / 2 ; \ k i 



1+2 



+k, 



T i + 2- T 
,h 2 /2 



Making use of these approximations in equation (2) gives the steady-state 
heat flow difference equation: 

[y T ] + (g +k )Mj =o 



. k 3 +k , 



referring to the general cell illustrated. 
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Exercises 



With this equation the standard technique of relaxation by residuals may 
be used to find the steady-state temperature distribution. 

This problem has been discussed in order to illustrate how the relaxation 
technique may be applied in cases of two added types of complexities: 

1. Cylindrical coordinate system 

2. Nonhomogeneous material 



1. Write a FORTRAN program to perform a matrix multiplication. The 
maximum size of the input matrices is 20 by 20 . The actual size of 
the matrix is L, m, n, and is to be read from data cards. 



Note: 



C ij = 



A ik • B kj 



k=l 



*2. Using the method outlined in the matrix inversion topic, find the coef- 
ficient inverse and solution for: 

x 1 + 2x 2 + 3x 3 = 3 

2x± + 3x 2 + 4x 3 = 5 

3x-£ + 4x 2 + 4x 3 = 5 

*3. The discussion of the problem on the iterative solution for homogeneous 
equations mentions the use of the iterative process for nonhomogeneous 
equations . 

Given the set of simultaneous equations 

AX = B 
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one such iterative procedure is to start with an approximation to the 
solution and solve for one of the xj using as values all other xj previously 
obtained. This procedure is repeated for each equation solving for a 
different x- in each. The entire procedure is then repeated until there 
is no significant change in the values of all Xy 

For the set of equations 

25x 1 + 2x 2 + x 3 = 69 

2x 1 + 10x o + x = 63 

JL 2 «5 

x l + x 2 + 4x 3 = 43 
the iterative form is 

x-l = (69 - 2x 2 - x 3 ) /25 

x 2 = (63 - 2x! - x 3 ) /10 

x 3 = (43 - x± - x 2 ) /4 

Iterate for the values of Xp x 2 and x~ which satisfy the equations. 

*4. Perform relaxation on the below grid which has boundary conditions: 

boundary y 

for the potential described by Laplace's equation 

} 2 V . a 2 v 



x y 



= o 



Make use of symmetry wherever possible. 



5 1 1 1 

4 

j 

2 

1 

12 3 4 
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SECTION VII: EMPIRICAL RELATIONSHIPS 



The previous chapters have discussed various ways of handling problems 
adequately described by theoretical equations. Some attention has been 
given to each of the general steps in computer solution of the problems, 
namely: 

1 . Selection of the appropriate relationships . 

2. Mathematical manipulations involving algebra or the calculus to 
arrive at a suitable evaluation equation. 

3. Choice of numerical technique when necessary. 

4. FORTRAN programming and actual computer evaluation. 

It must be recognized that this approach does not suffice in situations 
where step 1 is impossible or where the known relationships do not 
adequately describe the problem. 

The lack of theoretical knowledge may often be adjusted for by the use of 
empirical data taken by means of experiments and test runs. 

The discussion of the vehicle simulation in a previous section introduced 
the use of tabular data and an empirically derived equation together with 
theoretical equations in the description of the vehicle operation. 

The engineer is no stranger to the use of tables of data. Page upon page 
of collected data arises in almost every test of an engineered part or 
system. Conclusions are as often drawn from the consideration of tables 
of data as from theoretical data. It is important that this source of infor- 
mation be available, as is, for use with computer analysis. At the same 
time the possibility often exists for derivation of an empirical relationship 
to fit the data. This section will consider some of the standard proce- 
dures for: 

1. Use of tables of data as they exist. 

2. Derivation of empirical relationships from the data itself. 

The discussion will be broken into four areas, as shown in the chart below: 





FUNCTIONS OF A 
SINGLE VARIABLE 


FUNCTIONS OF 
MULTIPLE VARIABLES 


TABLE 
LOOKUP 


a. Table storage 

b. Table interpolation 


a. Table storage 

b. Interpolation logic 


DATA 
FITTING 


a. Fitting Criteria 

1. Selected points 

2. Least squares 


a. Linear regression 

b. Multiple linear regression 

c. Nonlinear Estimation 

d. Dimensional Analysis 
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Functions of a Single Variable 



Given the following table of data, how might the proper value of Y for a 
given value of X be obtained during a computer run? 



X 



X(I) 







0.0 


X(l) 


4.0 


X(2) 


8.0 


X(3) 


12.0 


X(4) 


16.0 


X(5) 


20.0 


X(6) 


24.0 


X(7) 


28.0 


X(8) 


32.0 


X(9) 


36.0 


X(10) 


40.0 


X(ll) 


44.0 


X(12) 



.913 
.930 
.941 
.946 
.948 
.950 
.951 
.948 
.944 
.938 
.928 
.914 



Y(I) 

Y(l) 

Y(2) 

Y(3) 

Y(4) 

Y(5) 

Y(6) 

Y(7) 

Y(8) 

Y(9) 

Y(10) 

Y(ll) 

Y(12) 



TABLE LOOKUP 

One procedure is to load the entire table into the storage of the computer, 
and then for a given value of X search the table for the corresponding value 
of Y. 

A FORTRAN program to accomplish the loading of the above table into 
storage and table lookup for several values of X is shown in Figure 14. 
Note that the search is accomplished with the use of an IF statement 
within a DO loop. For the case of the argument X being exactly equal to 
a table entry value of X, the corresponding value of Y is simply selected 
and printed. When the argument X falls between two table entries linear 
interpolation is performed. The graphic derivation of the linear inter- 
polation equation is: 



Y = Yi + 



Yi+l 



fri+i - yi) < x " x i) 

x i+l " x i 



yi 




>• x 
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This is the equation of the straight line joining points (xj, y^) and (x i+1 , yi +1 ). 
Evaluation of the right-hand side for a particular value of x gives the 
corresponding value of y. This equation is used in programming statement 4. 



1 c for comment 

Istatement ; 

* num.er 3 FORTRAN STATEMENT 

1 3 ]6 7 10 IS JO :5 30 35 40 45 50 55 60 65 70 7J 


c 
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> ■ i • ■ ■ i . . . . 
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1 
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50 TS. X . , ...... 












* 
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l))*(TX-X(j -i )/,(X/l)r X(l- 


^.)J 










0/ 70 .*..,....,...., . 












. . . S 


c^k. =r ! / /.> ...... 












8 
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G<i T4 1 ........ . 


. . . i .... . i . . 












AW. ........... . 


.. . . . i . . 












... 1 ... 1 .... 1 , 1 ... , 


... i .... i .... i .... i ... . 

























Figure 14. 

If the approximation of the function by a straight line in the interval 
(yj+1 ~ yj) is not of sufficient accuracy, the function may then be approx- 
imated by a parabola or higher-degree polynomial by such methods as the 
Lagrange interpolation formula. 

DATA FITTING 

If the problem involving use of this table were to be run many times on a 
computer, consideration may be given to finding an equation which will 
pass either through or within tolerance of all the points in the table. In 
most engineering problems it is sufficient to find an equation which passes 
within a specified tolerance of all the points in a table of data. 

Model Selection 

First a selection of the equation form to be fitted to the data must be made , 
A few of the possible selections are listed below. 



TYPE 



1 . Polynomial 



2. Logarithmic 

3 . Exponential 

4 . Power 

5. Fourier Series 



EQUATION 

From y = slq + a,x 

To y = slq + a-^x-L + a 2 x 2 + a 3 x 3 + a 4 x^ + a 5 x 5 

Generally restricted to this range. 

y = a + a x log x 

y = a^x 



y = a n x 



m 



y = a„ + 2_, (a cos nx + b sin nx) 
n = 1 
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The actual decision may be made on the basis of theory, some preliminary 
plotting, past experience, or by trial and error. 

Model Fitting Criteria 

Next a criterion for actual fitting of the equation form (finding values for 
the a^) must be selected. Some of the possible methods are listed below: 

Selected Points. As many sets of observation data as there are a^ to be 
determined are substituted into the selected equation, and the resulting 
system of equations is solved for the a i - While this method is very crude, 
it may be of value in situations in which available data is limited (see 
conical filter shell stress problem later in this section) . 

Harmonic Analysis . This widely used computer application is useful in 
fitting a Fourier series to a set of periodic data. 

Least Squares. This is the most widely used procedure for calculating 
the parameters a^ for the selected model. 

All of the models mentioned above (Fourier series excepted) may be fitted 
to a set of data by least squares. The principle will be explained with 
respect to the simple linear model 



y = a o 



+ a-, x 



Given a table of n sets of data, 
where y is assumed to be the above 
linear function of x, determine slq 
and a]_. 




The least squares criterion is to determine a Q and a-^ such that the sum 
of the squares of the vertical distance between the data points and the 
straight line is a minimum. Referring to the graph, this may be stated 
mathematically as: 



/ , i = mini 



minimum 



i = 1 
This is true only if 



> E s i 2 

i = 1 



3 a 







105 



and n 

> L € i 2 



i = 1 — =0 



3 a x 

The sum may be expressed in terms of the equation to be fitted and the 
original data points. 

€ 4 = Yi - (a + &! 34) 
Then: 



1^ u 

} L, € i 2 ^ Z^ [y. - (a + a 2 x.)] 

i = 1 i = 1 =0 



3a 5 a 



2 

5 Z_/ € i 2 ^ Z_y [yj - (a + a x 34) 1 

i = 1 = i = 1 _ 

'ia.i ~ > a 1 

Performing the differentiation and simplifying gives: 



na Q = 



11 ■" 

L-j y 4 - Z_, a i x i " 

i = 1 i = 1 

n n n 

T T 2 y 

L—t x i y. - a Z__/ Xj - a Z_^ x. = 
i = 1 i = 1 i=l 

These two linear, nonhomogeneous equations may be solved for a ft and a, . 
The parameters a.Q and a, , therefore, are computed in terms of sums and 
sums of cross products of the raw data. 



Functions of Multiple Variables 
TABLE LOOKUP 



Given a set of engineering data of the following form , from which x is to 
be computed for various sets of values of y and z: 
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x i yi 

x 2 Y2 z l 

x 3 Y3 

x 4 y_4 

x 5 ys 

x 6 YQ z 2 

x 7 Y7 

x _8 y_8 

x 9 Y9 

x io yio 

x n yn z 3 

A table lookup with interpolation procedure, more complex logically to be 
sure, may again be used. This procedure may be stated as: 

(a) For zi interpolate for x as a function of y alone. 

(b) Store the resultant value of x along with z^ in another table. 

(c) Repeat steps (a) and (b) for all values of z giving a complete table of 
x as a function of z alone . 

(d) Interpolate in this resultant table for final value of x to be used. 

DATA FITTING 

Quite often, particularly in such fields as chemical engineering, the table 
lookup possibility is impractical for one of two reasons: 

1. The problem demands a fitting equation for a meaningful solution. 

2. The data cannot be obtained in a form similar to that shown in the 
previous table. 

When this situation is true, curve fitting may again be used. 

For example, an equation in the description of the vehicle simulation 
model was presented as an empirical relationship. In this case R is 
assumed to be a function of three variables M, A and V. 

2 
R = aM + bAV 

The least squares criteria may again be used in conjunction with exper- 
imental data to determine values for a and b which best fit the equation to 
the data. 

To apply least squares relationships in this instance would require calling 
AV^ a new variable x 2 - AV , giving a linear relationship of the form: 

Y = ax 1 + bx 2 
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Statistics 



In general, data fitting of linear equations is called linear regression. 
If, as in this case, there are more than one independent variables, the 
procedure is called multiple linear regression. 

This chain of first assuming the form of a relationship and then using 
mathematical criteria to fit the relationship to experimental data can be 
thought of as a search for a useful "prediction equation." 

Note that in the discussion, the concept of "best-fitting" predictive 
equations to data has been used. The assumption has been made that a 
useful equation need not fit the complete data set exactly. This is an 
assumption based on statistical principles. 

An examination of statistics indicates further that predictions, or useful 
engineering conclusions, can be drawn from data without necessarily 
performing data fitting to arrive at an equation. The most common sta- 
tistical topics will next be discussed. 



The general topic of statistics has been perhaps the most neglected of all 
engineering tools for industrial work. A major reason for the neglect has 
been the vast amount of calculation necessary to obtain all the possible 
statistical information from a set of data. This "excuse" is no longer 
valid with computers available to do the necessary calculation. Futher- 
more, the most important statistical procedures have already been 
programmed for many computers. To achieve rather complete statistical 
reduction of information, the engineer need only prepare the data in a form 
acceptable to one of the available programs . 

It should not be necessary to point out that the design of the experiment or 
test should depend upon the statistical procedures to be applied to the 
resulting data. Before any data is taken, the one or more engineering 
hypotheses which are being tested should be clearly defined. The sta- 
tistical procedures available to resolve the test should be determined. 
Once these are decided upon, the design of the experiment to supply the 
necessary data can be most profitably begun. 

Techniques, which will not be discussed here, have been programmed to 
assist in the determination of the number and order of data required in an 
experimental design. 

Statistics is a form of shorthand capable of reducing volumes of data to 
lines or pages. These statistics are descriptive in the sense of making 
positive statements about the meaning of the collected data. Furthermore, 
the statistics indicate the confidence which can be placed in any statement 
which is made, the accuracy to be expected, the range of application of 
the statement, and the probability of the statement being true. Decisions 
can be based on more complete knowledge. 

Levels of Statistical Evaluation 

There exist a great number of statistical methods which may be applied to 
any given set of data. At the same time data may arise in an infinite number 
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of forms and quantities. Interest may center on the measurement of one 
variable, as in the expected life of a brake shoe under test, or on the 
intereffects of several variables, as in the study of factors determining 
sales of a product. The more common statistical methods and the type of 
results desired fall naturally into a pattern which requires a progression 
from simple to complex calculation procedures. Four levels in this 
progression may be identified as: 

1. Probability analysis of a random variable, which is useful in quality 
control, testing of vendor products, and field performance of products. 
No cause and effect is considered. 

2. Analysis of variance to determine the significance of differences be- 
tween classed or grouped data, such as the differences caused by variation 
in the process for preparing a product. This represents a test for the 
existence of cause and effect. 

3. Correlation analysis , which gives a measure of the linear relationship, 
dependence or association between two variables. As such, it represents 
an attempt to place a measure on the cause -and -effect relationship. 

4. Regression analysis is a computational method for determining para- 
meters in an assumed equation form which expresses the dependence of 
one variable (the dependent variable) on one or more other variables 
(independent variables) when data on all variables is available. This, then, 
is a method for defining the cause-and-effect relationship to the extent 
that useful predictions can be made for the behavior of the dependent 
variable. 

Probabilities 

For many problems in which questions of probabilities arise, the behavior 
of individual events in the system is known beforehand or a priori. For 
example, the tossing of coins, the throwing of dice and the appearance of 
cards in a dealt hand are ruled by laws of probability which are well de- 
fined . 

The laws for dealing with a priori probabilities will not be considered here. 
In engineering work the a priori probabilities are rarely known, but the 
probabilities for an event can be determined on the basis of data obtained 
by experiment or testing— for example, the expected life of brake shoes. 
This is a posteriori probability and represents an estimate (a useful esti- 
mate) of the true probability. The estimate can be improved by increasing 
the accuracy and number of trials made in the experiment or test. 

If every possible trial is made (that is, if every brake shoe to be used is 
tested) the "population" has been tested. However, it is generally practi- 
cal to test only a "sample" of the population (that is, many tests are 
destructive). From this testing of population, or sample, more than just 
the probabilities of a single event is desired. What is desired is an 
inventory of all possible values for the event and the probabilities of the 
event taking on each value. This inventory is called a probability 
distribution. 
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Discrete distributions are used to describe probabilities for which events 
can take on only discrete values (that is, the number of heads in tossing 
coins). The binomial distribution and the Poisson distribution are ade- 
quate to describe most naturally occurring discrete distributions. These 
will be ignored here since the majority of engineering problems involve 
continuous distributions (for example, the life of a brake shoe which might 
take on any time value dependent on testing procedure and accuracy of time 
recorded). 

This discussion will be restricted to the normal distribution for continuous 
distributions of probabilities although the Chi-Square distribution, the 
F-distribution and the t-distribution should certainly be considered in any 
extensive statistical investigation. Statistical tests can be made to de- 
termine the suitability of application of the normal distribution to a 
particular set of data. 

To illustrate the use of probabilities and the normal probability distribution, 
consider the table below showing hypothetical results of the life tests for 
a particular brake shoe . 



Life in hours Xj 


Frequency Fj 


Normalized 




(Midpoint of Interval) 


(Number failed in Interval) 


Frequency fj 


F-x- 
11 


10 


50 


.050 


500 


20 


75 


.075 


1500 


30 


100 


.100 


3000 


40 


200 


.200 


8000 


50 


250 


.250 


12500 


60 


150 


.150 


9000 


70 


75 


.075 


5250 


80 


75 


.075 


6000 


90 


25 


.025 


2250 


Sum 


1000 


1.0 


49000 
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A common method for displaying this information is to construct bar 
graphs to show the frequency distribution and another to show the accumu- 
lative distribution. 



1.0 
.9 

.8 
.7 
.6 



Graph A 1 

Frequency 
Distribution 
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Graph A 2 
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k 



10 20 30 40 50 60 70 80 90 100 
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The assumption can be made that for successive decreases in the interval 
size and increases in the sample size, the above graphs will approach the 
continuous curves as shown below: 



Graph B 1 



Frequency 
Distribution 



Graph B 2 



Accumulative 

Frequency 

Distribution 




90 100 




10 20 30 40 50 60 70 
Life in hours Xj 



90 100 



10 20 30 40 50 60 70 
Life in hours Xj 

These graphs, or the listed data, may be consulted for the answers to 
such questions as: 

9 What is the probability of an individual brake shoe having a life greater 
than or less than a given value of x? Quality control may dictate an ac- 
ceptable probability of the life being less than a certain number of hours. 
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•What is the probability for the life of an individual brake shoe lying between 
two values of x? Again the variance is of importance in quality control. 

A demonstration of the shorthand quality of statistics is the fact that the 
information displayed above can be captured in two statistics, assuming 
a normal distribution: the mean and the standard deviation. 

The mean of the sample is 



N 



N 



where N = number of observations 
for ungrouped data 



or 



N /M 



where M 



number of intervals 
for grouped data 



The mean x is an estimate for the mean of the population (ji). In computer 
handling of data, the grouping of data — a device for simplification of 
calculation— is seldom used. From this point ungrouped data only will 
be considered. 



The standard deviation of the sample is 

-J4 



S = 



N 



x) 2 



(N-l) 



and s is an estimate of a, the standard deviation of the population, 
called the variance. For a normal distribution the probability of an 
individual reading lying within X + cr is approximately .68. 



s 2 is 



With the use of X and s, the graph b^ can be constructed from the use of 
the normal distribution function. 



f(x) 



= svW 



-1/2 ( x - xV 
e s~ 
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This function has been tabulated for the "standard" normal distribution where 
the transformation 

+ _ x - x. 



gives the distribution function 

2 

-t 



1 2 

f (t) = G 
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Graph b2 can be reconstructed from the integral function 



(J) (x) = 1 - / f (x) dx 



to give the probability of an individual reading being greater than X. 
Again, this function has been tabulated for the variable 



_ x - x 



or 



t = 



<j) (t) = 1 - / f (t) dt 



For hand calculation a table of the general form 

t f(t) (|)(t) 

-2 .05399 .0228 

-1 .24197 .1587 

.31894 .5000 

1 .24197 .6587 

2 .05399 .9772 



can be used, along with the transformation equation t= - x " - to answer 
frequency and probability questions. 

For example, with a mean X = 5 and a standard deviation s = 2, the proba- 
bility of an individual reading lying between 1 and 9 can be calculated as 
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l l~ 2 ~ " 2 1 2 2~ ~ 2 




P (l< x^9) = / f < x > dt - / f (t)dt = 4> (ta) - <p (t x ) 



= .9772 - .0228 
= . 9544 



However, in computer handling of probability problems it is simpler to 
calculate values for f (t) and (p (t) than to store the large table in the 
computer. The integral evaluation for cf> (t) is difficult because an exact 
integration cannot be performed. The Hastings approximation * for 



y 2 

e~ Y dy 



o 

is often used with appropriate transformations to evaluate the integral in 
a computer program . 

Confidence in Mean and Sample Size 

Questions concerning the confidence which may be placed in the calculated 
mean and in the chosen sample size can be answered by consideration of 
the statistic 



x -m 



which has the Students t-distribution. The Students t-distribution 
approaches the normal distribution as a sample size increases and in 
this discussion a normal distribution will be assumed. 

To compute the "range on the population mean for a given confidence", 
t_ is chosen according to the confidence (in the table above, for t = 2, 
the confidence that the population mean will be such that the calculated 
t c = 2 is approximately .97). 



* Approximations for Digital Computers , by Cecil Hastings, Jr. , Princeton University 
Press, 1955. 
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If 

U = X + A 
then, using (1) above 



A = t c S (2) 



and the statement can be made with a 97% probability of being correct 
that from the sample data, the mean of the population lies within 

jU =X + A 

Equation (1) may also be used to determine whether a sufficient sample 
size has been chosen to give a set confidence in the population mean lying 
within an acceptable percentage variation, K, of the calculated sample 
mean. 



If 

\X = X ± KX for 

then from (2) 



KX= t c S 



£ KX 



yNc" 



and 



2 2 
N c = C _p (3) 

(Kxr 

where N is to be compared with the actual sample size N. If N c > N, a 
larger sample will be required for the set confidence. If N c < N, a 
sufficient sample size has been used. It should be noted that because of 
approximations made, and arbitrary choice of initial sample size, the 
value of N only indicates the direction in which the sample size should 
be changed— not the actual size change required. Several iterations (new 
sample size - calculations of N ) might be necessary to determine a best 
value for N for a given class of data. 

ANALYSIS OF VARIANCE 

Si 2 2 

The statistic F = gSr » where Sjl and S2 are variances of samples from 

populations with normal distributions whose true variances are equal 

(i. e. , O-]^ = O 2 , which again can be tested), has a distribution of the shape 

shown below. 
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g(F) 




The F test is used to test whether there is a significant difference in two 
sample variances. 

The test consists of setting a confidence level (as in the discussion above 
on the mean, using the normal distribution) in terms of percent of area 
(shaded area = ~< and total area = 1) which may lie to right of F,^ in 
distribution curve. This determines the Fc: and means that any given F 
greater than F^ has only a probability <* of occurring due to random 
chance alone. The confidence in F being less than Fo< is 1 -•< . 

The distribution function g (F) is a complex multivariate one (dependent 
on F, degree of freedom of Si and degree of freedom of S2 ) and for this 
reason tables are generally tabulated only for °< = .01 and «*« = .05. For 
the same reason g(F) is not as often calculated as part of a computer 
program as is (|)(t) for the normal distribution. 

2 2 

Next the sample variances S-^ and S 2 are calculated and F observed is 

computed. If F>Fo< there is confidence of 1 -«=-< that a significant differ- 
ence in the sample variances exists. 

o 
For analysis of variance, S^ is a measure of the variation in test data 

caused by a difference in the process or treatment, while S2 is a 

measure of the purely random variation of the test data. If the resultant 

F is greater than the preset F^ , then an effect on the data results can be 

attributed to the variation in the treatment. 

2 
The assistance a computer can afford in calculating the variances S ^ and 

S 2 and the resultant F can be seen from the following data format and 

required calculations for an analysis of variance. 
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Analysis of Variance Data 
Replications 



Treatment 


1 


2 


3 


. 


n 


1 


x i,i 


X l,2 


• 


. 


x i „ 
1 > n p 


2 


X 2,l 


X 2,2 


• 




• 


3 


X 3,l 




• 


• Xp.i • • 


* 

X 
P,n p 


k 


^.i 




• 


• 


\n p 



* All n^ are not necessarily equal - there may be a difference in the 
number of replications for different treatments. 



F observed =■ 



Calculations 

2 9 

(S^ has k -1 degrees of freedom and S2 has 

k 

S (n n -l) degrees of freedom) 
p=l p 



where 



k n p 



p=l i=l P' x / p=l p 



k n p 

TSS = E, J x . - C 
P=li=l P.i 



2/ n 2 2/ n k , 

A + ia^i v n 2 + - - - - is x k,i / n k - C 



ni 2. n 2 2 

Gss= & x i,i>A + ia x 2,i ) / n 2 + 



and S^ = GSS/(k-l) 



3 2 = (TSS - GSS)/r (np-1) 
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A sample analysis of variance calculation performed on a computer is 
shown below: The observed F ratio of 15.22, when compared with tabu- 
lated F (degree of freedom for S^ =8, degree of freedom for S2 =65, 
which is closest entry to 63) where Fo< = 2. 08 for «x = . 05 and F«< = 2. 79 
for <=<= .01, shows the variation in data due to the variation in treatment 
to be significant at both the 5% level and the more stringent 1% level. 
The treatment has a cause-and-effect relationship with the variable for 
which the data was recorded. 



TREAT- 


NO. OF 


REPLICATIONS 












MENT NO. 


REPS 
















1 


8 


6.60 


5.11 


7. 1U 


U.89 


6.81 


U.93 


U.33 


2 


8 


2.67 


2.19 


2.36 


0.31 


1.U8 


0.68 


1.19 


3 


8 


5.90 


U.58 


2.08 


3.07 


3.12 


3.73 


U.00 


h 


8 


5.80 


5.92 


U.88 


5.32 


5.36 


U.6U 


5.68 


5 


8 


5.00 


5.58 


6.08 


U.U1 


5.37 


5.89 


6.3U 


6 


8 


2.95 


U.58 


2.89 


3.93 


1.67 


U.61 


3.28 


7 


8 


5.86 


6.1U 


5.93 


U.23 


U.59 


5.18 


5.60 



U.80 6.U7 5.63 U.68 



5.U2 6.93 3.62 5. 1U 



TREAT- 


SUM 


SUM 




MEAN 




MENT NO. 


Y 


Y2 








1 


U5.68 


268. U5 




5.71 




2 


13.78 


30.07 




1.72 




3 


28.68 


11U.02 




3.58 




U 


U0.06 


209.38 




5.01 




5 


UU.81 


253.98 




5.60 




6 


26.66 


95. 8U 




3.33 




7 


U2.6U 


230. U7 




5.33 




8 


U2.56 


231.00 




5.32 




9 


U3.72 


2U9.78 




5.U6 




ANALYSIS OF VARIANCE 








SOURCE OF 


DEGREES OF 


SUM 


OF 


MEAN 


VARIATION 


FREEDOM 


SQUARES 


SQUARE 


AMONG 




8 


120 


.8U 


15.11 


TREATMENTS 












WITHIN 




63 


62 


• 5U 


0.99 


TREATMENTS 












TOTAL 




71 


183 


.38 





The above discussion has been related to data with a one-way classification 
according to treatment. Data may also be classified according to two or 
more different treatment variations. The experimental design and analy- 
sis calculation will vary accordingly. 

CORRELATION AND REGRESSION ANALYSIS 

A brief discussion of the computer applications in the calculation of corre- 
lations and regression equations is simplified by the use of the matrix 
notation introduced in a previous section. 

This discussion will concern the information to be gained from test or 
experimental data which is available in the following form. 



Observation 


Dependent 


Independent Variables 






Variable 










Y 


x i 


x 2 x 3 


X P 


1 
2 
3 


^1 
^2 
^3 


X 2,l 


X l,2 X l,3 ' 
X 2,2 X 2,3 


X l,p 


N 


% 


Xnl 




X n,p 
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Regression analysis consists of using the above data to determine according 
to the least squares criterion the values of bj which best fit an equation of 
the form 

y = b 1 x 1 +b 2 x 2 - b p x p 

to the data. The objective is to obtain a useful prediction equation. 

The value of P must be less than N for a correct analysis. If a constant 
term is desired in the equation— that is 



y = Vi + V2 + b p-l X p-l + b p 

a column of ones would replace the Xp column in the above data. 

The least squares analysis for this case can be most simply described 
using the matrix handling rules of an earlier section. 

Let 

77 = observed column matrix of observed y values. 

B = parameters bi to be determined, a column matrix. 

D = rectangular matrix of observed values for the independent variables, 

and set 

rj = DBor (1) 



^1 = X l,l b p + X l,*2 + x l,iJ>p 

?2 = X 2,l b l + X 2,£ 2 + X 2,p b p 



77 = X Jb, +X J3 + X b 

'n n, 1 1 n, z 2 n,p 



n,p p 

where n > p. These equations may not be solved for B since D has more 
rows than columns. However, multiplication of both sides of the above 
equation by the transpose D 7 of the D matrix gives 

d 1 n = d'db < 2) 

where D D is a square matrix. These equations are equivalent to the 
summation equations which resulted in the least squares analysis, called 
the normal equations. 

The solution of these equations for B is found by first obtaining the inverse 
of the D D matrix and multiplying both sides of (2) by the inverse (D D) -1 
to give 

(D 7 D)" 1 d't? = (D* D)" 1 (d' D) B 
or B = (D 7 D)" 1 D ; 77 (3) 
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The several matrix manipulations indicated in this solution for the values 
of b. imply so much computation that for any problem involving four or 
more independent variables a solution without a computer is virtually 
impossible. 

STATISTICAL EVALUATION OF A REGRESSION ANALYSIS 

Once the evaluation of the bj by equation (3) has been completed, several 
statistics become available to judge the value of the analysis for pre- 
diction purposes. 

• "Goodness of fit" - the standard error of the estimate 

Let Yj = the predicted value for y using the original data 



yi = b l x i,l + b 2 x i,2 + b p x i,p 

m= the observed value for y 

ei = Yi - Vi 



Then the standard error of the estimate is given by 



2 N 2 / 
S e = ^ e i )/(N-p) (1) 



a statistic analogous to the variance for a single random variable (that is, 
68% of predicted values should be in error by less than + Se). 

• Simple correlation 

One method for computing the simple correlation coefficients between 
each possible combination of two variables at a time is by 

r i,J = V i.j/ s i S J 
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where 

V i»j = ~W < 2x i x j _ x i 2 x j ) 

s i = (^ i) ' and y can be considered one of the variables 
X[ or Xj. 

The interpretation of the simple correlation coefficient between two 
variables is that the square of rj j is the percentage of the variance 
of xj. that is accounted for by its relationship with Xj. This applies only 
if a linear relationship can be assumed, and it ignores the possibility of 
intercorrelations with other variables. 

• Partial correlation coefficients 

Let aj j be the elements of the inverse of the simple correlation coefficient 
matrix. Then 

x • = ~ a y> j 



tj a y>y * a JJ 

are the partial correlation coefficients. The interpretation of this statistic 
is that it gives the true correlation between each pair of two variables (one 
must be the dependent variable) out of the total investigated after the effects 
of the remaining variables have been taken into consideration. Intercorre- 
lations can make this vary a considerable degree from the simple corre- 
lation coefficient for the same two variables. 

• Multiple correlation coefficient 

A measure of the total variation of the dependent variable y that has been 
accounted for by the regression analysis, analogous to the simple corre- 
lation coefficient for two variables only, is called the multiple correlation 
coefficient and is given by 



R 



■^ 



l y>y 



This gives an excellent measure of the success of the regression analysis 
approaching 1 as the regression equation improves in fitting the data. 

• Alternate calculations of the standard error of the estimate 

The calculation of S e in (1) above implied a considerable amount of ad- 
ditional calculation. Two additional measures, more simply calculated 
once <ajii is available are: 
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Biased Standard Error of the Estimate 



- s y 

^'e ~ rs , and the 

V a y,y 



Unbiased Standard Error of the Estimate 



VN-p 



P+l 



Significance test for independent variables 



The F test is also used to test the hypothesis that a regression coefficient 
bj is zero (one of the variables X^ can be ignored). 



Hypothesis b^ = 

Test: 



b.. 2 



F observed = F^ n _^ = — l 



,a i,i 



where the subscripts on F indicate one degree of freedom for the numerator 
and N-k degrees of freedom for the denominator. 

If F observed < F^ then assume bj = or that Xj can be removed from the 
analysis. Fo< can be set for either the 1% or. the 5% level of significance. 

Study of Anodized Aluminum Panel Corrosion 

This example of the use of both analysis of variance and regression 
analysis will illustrate the application of the above-described techniques . 

A process engineering group is asked to evaluate 60 different processes 
or treatments for anodizing aluminum panels to be used in the body of an 
automobile. One lab test A has been the standard for evaluating the cor- 
rosion resistance of anodized panels, but now a new lab test B is proposed 
as being less expensive to perform and perhaps more accurate. It is 
assumed that the environment has an effect on the corrosion and that the 
thickness of the panel has an effect, but the relative importance of each 
is unknown. The decision is made to conduct a statistical study of data to 
be collected on the corrosion of aluminum panels of various treatments 
which have been subjected to a variety of environmental conditions. 
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The objectives of the experiment are to: 

1. Measure the relative influence of environment on corrosion. 

2. Measure the relative influence of panel thickness on corrosion. 

3. Determine whether a significant difference in corrosion can be traced 
to the variation in treatments. 

4. Rank the treatments if a significant difference exists. 

5. Evaluate the relative merits of lab tests A and B as predictors of 
corrosion resistance. 

6. Determine whether the corrosion resistance of a particular panel can 
be predicted in terms of one or both of the lab tests, the thickness of 
the panel and the environment. 

Consideration is given to the statistical tools available to analyze the data. 
It is decided that ten panels shall be placed on each of 100 cars which are 
used in widely varying environmental conditions . Nine of the ten panels 
are selected at random for each car from the 60 different panel types, so 
that 15 panels of each treatment are used. The tenth panel is a control 
panel for the establishment of an environmental rating. The control panel 
on each car has the identical treatment and the same thickness, so that the 
only variation (as close as it can be controlled) in its corrosion resistance 
will be due to the variation of environment of cars . 

After the panels have been on the cars for four months they are removed 
and a relative measure of their corrosion is made on the basis of accepted 
industry tests. Since the scale of corrosion rating (1. to 10. 0) has been 
chosen on the basis of past experience, it is expected that the average 
rating (called R) will be 5.0, with a standard deviation of about 2. 0, so that 
for a normal distribution panels are expected with rating from 1 to 10.0. 
Although difficult, the relative accuracy of ratings is kept to about 2% since 
it is recognized that the statistical analysis can give results of no greater 
accuracy than the data used. 
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The result of this data collection is a set of data with the form shown in 
the table below. 

Anodized Aluminum Corrosion Data 
Control Data 



Car 
3 4 



100 



Resistance 
R 



3.0 7.6 5.5 4.3 



9.0 



Treatment Data 



Treatment 
Resistance 

1 R 
Thickness 

T 
Environment 

E 
Lab Test 

A 
Lab Test 

B 

2 R 

T 

E 

A 

B 

R 
t 



60 



i 
t 

t 

R 



Replication 
1 2 3 4 

6.8 5.9 

.015 .014 

3.6 7.8 

35 31 

39 34 



15 
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Using this data with computer programs: 

1. The control data for resistance R was given the probability analysis 
for a random variable and an environment rating was set up (ranging 
from 1.0 to 10. 0) so that 10% of the cars fell into each range of 1. 0. 
The data on environment for each car was then filled in on the general 
treatment data sheet. 

2. Ignoring the data on thickness, environment (they are assumed to fall 
at random) and the lab tests A and B, an analysis of variance was made 
for significance of difference in the variable R from one treatment to 
another. A high level of significance was found. 

3. The treatments were ranked according to their means with a secondary 
consideration given to the standard deviations of each treatment. 

4. A regression analysis was made assuming the equation form 



R = b 1 t + b 2 E + b 3 A + b 4 B 

The relative influences of the thickness and of the environment were de- 
termined by an examination of the simple correlation coefficients, the 
partial correlation coefficients, and the F tests for the significance of 
factors . 

The same tests were applied to determine the relative merits of lab tests 
A and B in predicting the resistance of R. 

The multiple correlation coefficient gave a measure of how completely the 
variance of R was accounted for by the prediction variables chosen. 

Finally, the usefulness and accuracy of the prediction equation was de- 
termined by the calculation of the standard error of the estimate. 

The preceding material has been designed only to whet the appetite of the 
reader for further examination of the potential values in statistical evaluation 
of experimental data now that computers exist to remove the burden of 
calculation. The growing use of statistics for practical industrial engi- 
neering work is reflected in the literally hundreds of computer programs 
now available for data reduction. 



Dimensional Analysis 



This technique for model selection of relationships between physical 
variables has wide application to engineering problems. With electronic 
computation the application of dimensional analysis can be greatly facil- 
itated. 

As usually written, each physical variable of an engineering equation is 
expressed in some combination of fundamental units . These units — time , 
length, mass, temperature, etc. — are the dimensions. 
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Dimensional analysis is a systematic study of the dimensions of an equa- 
tion containing several variables in order to: 

1. Discover any variables whose dimensions are incompatible. 

2 . Ascertain the nature (but not the true value) of the relationship between 
the variables. 

3. Assemble the variables in dimensionless groups so that the relation- 
ship between these groups can be expressed algebraically. 

As an engineering application of dimensional analysis which will make the 
procedure clear, consider a problem of conical filter shell stress. A 
filter shell must be designed to withstand large stresses applied through 
a tightened nut at the top. This force is needed to overcome the force due 
to oil under pressure inside the shell and maintain a firm seal at the base 
gasket. 



) Reinforcing 
Ring 




w 



Shell 



IHJ~| 



Base 
Gasket 



An equation is desired which will predict the applied stress, W, for which 
the shell will fail. For this structure a complete theoretical analysis is 
too involved. Therefore, an empirical relationship will be determined. 
It is assumed that variables involved are: 

S = Ultimate shear stress of the material 

t = Thickness of the shell 

D = Mean diameter of shell 

= Cone angle with the vertical 

d = Diameter of reinforcing ring 

W= Applied stress 
Dimensional analysis now consists of finding dimensionless groups of 
variables so that a dimensionally correct equation may be written. An 
approach here is to let 



7T ± = S a t b D c d 
7T o = S e t f D g 



(1) 



7T 3 = S h t 1 DJ W 
The dimensions of each of the variables to be inserted are: 

S = MT' 2 L" 1 

t = L \ M for mass 

D = l "S L for length 

d = L ( T for time 

0=1 

W = MT" 2 L" 1 
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For it i to be dimensionless, the sum of the exponents for each funda- 
mental dimension must equal zero. 

For M, a = 

For T, a = 

For L, -la + lb + lc + 1 = or b = -(1+c) 

For this particular case the resulting three equations have an infinite 
number of solutions. For simplicity consider: 

c = 
Then 

b = -1 

and 7T x = S t _1 D° d 

or7r= (A.) 
1 t 

Similarly, tt 9 = (-H ) e 

& t 

(Since 9 is dimensionless, this tt factor may be rewritten as two dimen- 
sionless groups.) And 

TT =(^) 

3 W 

At this point the variables are arranged in dimensionless groups and a 
relationship between the groups may be assumed. 

A reasonable form to assume is: 

7T 3 = K(7r 2 ) a (Tr^b 
or 

st 2 = W_D_\ a /A\b e c (2) 



t -(f) (ty 



where 9 has been used as a separate it factor because it is dimension- 
less. 



Rewriting (2) gives 

-f-(fr (f) be 

where K, a, b and c are constants to be determined 



KW /_D\a /_d_\ b c 

t 2 it/ ItJ (3) 



One method for calculating suitable values for the constants K, a, b and 
c is that of selected points covered earlier in this section. 

(a) Take the log of each side of (1) 

log S = log (~2-) + a log (— ) + b log (-^-) + c log 9 (4) 
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Probabilities 



(b) Make four shells of different dimensions with known ultimate stress, S. 

(c) Crush shells to obtain W. 

(d) Substitute all the values into (4) to obtain four equations in four un- 
knowns (K, a, b, c) and solve these equations simultaneously. 

(e) Apply constants to (3) and verify by testing other shells. 

For improved accuracy at the expense of shells, the following procedure 
may be used (this is equivalent to the least squares done graphically): 

(a) Holding S, t, D and d constant, crush many shells of varying . 
For this case (4) becomes 

log W = c log + B 
where B is a constant. 

(b) Plot W versus on log-log paper. The slope of the straight line is 
then c. 

(c) Repeat a and b for other variables . 

Note the possibilities for computer usage on larger problems of this kind. 
Simultaneous equations (to be discussed later) are handled twice. Further- 
more, the procedure for improving the accuracy is a simple curve fitting 
problem. In fact, the whole procedure for dimensional analysis, since 
it is well defined, can be programmed for a computer. 



The probability of a variable having a certain value is a question which 
often arises in engineering. If historical data on the variable is available, 
it is often possible to do the necessary predicting. This problem illus- 
trates a method for handling normally distributed variables . 

Given a table of observed values for a particular variable whose frequency 
distribution is assumed to be normal: 

(a) Find the probability of a reading being less than any given value of 
x, and 

(b) Find the probability of a particular reading lying between any two 
values of x. 
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The general distribution function is: 



l/x-m\ 2 



f(x) = 



2\ <r 



2w 



where m = mean 

a = standard deviation 



The standard normal curve is given by: 

-t 2 /2 
e 

f (t) = ;= 

V2 7T 

Therefore a transformation of variables will allow us to make use of 
standard normal distribution function. This transformation is 



t = 



x-m 



The probability of a reading falling below a given value of x is given by: 

/ f(x)dx=/f (t)dt 
-co 



-00 



x-m 



where t =■ 



The probability of a reading lying between x^ and X2 is: 



r x 2 r 2 

J f(x)dx = J f(t)dt 



x n -m 
where t, = — 



t = 2 

l 2 — 



Xo-m 
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Exercises 



The problem, then, is one of evaluating the integral 



b -t 2 



hi e 2 dt 



V~2tt 
a 



A table of values of this integral for various values of t (normal areas and 
ordinates) is available in any statistics text. In computing, problems of 
this type might be tackled in two ways: 

1 . Store the table in memory . 

2. Use numerical approximation for evaluating the integral. 

An approximation by Hastings* might be used with the latter approach. It 
is 



2 
o 



x 

.2 



~ 2 r -t" 
<P ( x ) = y^r J e dt for o < x < oo 



1 

where (J) (x) = 1 - -.4 

[l+a^x+agX +a^x +a^x J 



and 



2l = .278393 a_ = .000972 

1 o 

a n = .230389 a. = .078108 

2 4 



1. Write a FORTRAN-language program to perform table lookup with 
linear interpolation for a tabulated function of two variables as 
described under "Functions of Multiple Variables." 

2. Evaluate the coefficients for the linear least squares fit 

y = a + a x x 



for the following table of data 



♦Approximations for Digital Computers by Cecil Hastings, Jr. , 
Princeton University Press, 1955. 
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14.0 64. 

13.9 54. 

13.5 55. 
13.4 56. 

12.6 46. 
12.6 51. 
11.8 42. 

11.4 47. 
11.3 48. 

10.5 35. 
10.3 44. 
10.3 29. 

9.3 31. 

7.6 14. 

*3. The distribution in measured diameter of the dirt particles normally 
deposited in a carburetor sediment bowl is a normal one. 

For a mean of .0025 centimeters and a standard deviation of .0002, find 
the probability of a dirt particle being greater than . 0027 centimeters in 
diameter. 

Use any statistics text for necessary tables. 
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SECTION VHI: OPERATIONS RESEARCH 



The reason for including this as a topic of engineering analysis for com- 
puters is that often the technical personnel of a company will be turned to 
for answers concerning the mathematical aspects of Operations Research. 
Books have been written on each of the subjects to be mentioned. The 
intention here is merely to present examples in such a way as to make 
possible the recognition of problem types as they occur. 

Operations Research can be thought of as the application of scientific 
methods, techniques or tools in the management decision area. It attempts 
to present a clear picture of the system under study and at the same time 
predict the consequences of alternate sets of actions — in other words, 
show the best decision. Methods for attaining either or both of these 
goals require knowledge of the system. Managed systems are, in general, 
complex, which means that vast amounts of data are required to describe 
these systems . Manipulation of this data to attain the stated goals gen- 
erally requires a computer. 

OR seems to be at the "clear picture" stage at the present time. The more 
complex problem of "prediction" has been hindered by the fact that 

1. The mathematical and logical methods employed by OR are unfamiliar 
to those most intimately concerned with the problems (the executives) . 

2. The historical data describing the system does not exist in proper or 
even intelligible form. 

However, consideration of the prediction level of OR while working ^n tne 
clear-picture level has the following benefits: 

1. Creating a clear picture of a system results in the organization of 
data into intelligible form. This data can later be used at the pre- 
diction level of analysis . 

2. Consideration, even without use, of more advanced techniques often 
gives a better understanding of, or "feeling for, " a system. 

Almost any procedure for presenting data to an executive can be thought 
of as part of the clear-picture operation. However, the computer allows 
concise data handling in many new areas at reasonable costs . Some of 
these are: 

1 . Graphic presentation of parts failures in the field , showing effect of 
engineering changes on failures, failures per unit of operation, fail- 
ures as a function of production group, etc. , on a day-to-day basis. 

2. Presentation of historical data concerning use of spare parts in the 
past so as to gain some insight into what an "all-time production" of 
a spare part should be. 
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3. Creation of a weekly report to a supervisor, showing jobs assigned to 
his shop, estimated time for each job, group within shop to which it was 
assigned, actual time for completed jobs, and year-to-date figures for 
number' of jobs assigned to each group with totals for estimated and actual 
job times. This would allow analysis of each group's performance. 

4. Machine loading calculations are now commonly being performed on 
computers to give the total load in hours for each facility in a plant 
for a given production schedule. This can quickly indicate need for 
additional facilities or, conversely, excessive idle time. 

5. Dramatic improvements have taken place in the area of determining 
parts requirements for finished assemblies through the use of com- 
puting equipments. This allows rapid response to changes in orders 
for finished assemblies. 

What are some of the prediction techniques presently being employed? 

1. Linear Programming: This technique may be used on any system 
which may be represented by a set of linear equations . 

There may be more variables than equations and therefore an infinite 
number of solutions to the set of equations. The object of linear pro- 
gramming is to find a particular solution which will optimize (make 
maximum or minimum) another linear function (usually a profit or 
cost function). This technique has been very successful in the blend- 
ing of petroleum products and the mixing of grain feeds . It is not un- 
usual for the mathematical model of a gasoline refinery to be as large 
as 100 equations in 500 variables. The actual mathematical techniques 
employed for linear programming will not be discussed here. However, 
it is worth noting that the techniques have been refined to the point 
where their application can be rather mechanical. That is to say, a 
knowledge of the mathematics is not required for successful use of 
linear programming. 

An important special case of linear programming is the transporta- 
tion problem. It may be used to optimize any system for which a com- 
modity is available at several sources and is required at several 
destinations, and for which the transportation rate for each possible 
source-to-destination combination is known. The most economical 
shipping assignment will be determined. 

2. Forecasting: Historical data is analyzed statistically to enable simul- 
taneous forecasting of many interrelated factors of management inter- 
est in terms of other known variables that affect them. Although, in 
general, the statistical techniques are not new, the computer allows 
consideration of many more variables and, therefore, the handling of 
complex management systems. 

3. Critical Path Technique for Scheduling: This technique has developed 
rapidly into one of the outstanding examples of successful application 
of Operations Research techniques. The technique is an answer to 
many of the problems associated with the scheduling and control of 
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projects— for example, the engineering design of a system, the con- 
struction of a large building, or the setup of a manufacturing process 
for a new product. 

The technique rests on the principle of separating the planning and 
scheduling functions in the control of a project. Planning for the 
critical path technique of scheduling consists of determining all the 
distinct jobs necessary to complete the project and determining the 
sequence in which these jobs must be performed. The major new 
device employed here is the construction of an arrow diagram in place 
of the traditional Gantt or bar charts. For example, the arrow dia- 
gram shown below might represent the jobs (A to G) required in a 
design of a new carburetor (significant scheduling problems have arrow 
diagrams many times the size of the one shown). The numbers in 
brackets indicate the length of time planned for completion of that 
particular job. 




Start i \ *>l I 'wwwwwvswwwwwtmiV 1 Finish 



The major advantage of this type of diagram is that it indicates all the 
interrelationships of the jobs in terms of jobs preceding, jobs following, 
and jobs proceeding simultaneously with each given job. 

The scheduling phase consists of determining the critical path (indicated 
by wavy arrows), calculating leeway or float for the noncritical jobs 
and, finally, determining on the basis of cost considerations and man- 
power availability, a means for reducing the length of the critical path. 
It is in the scheduling phase that the use of a computer becomes 
imperative for projects of a significant size. 

There exist a great number of computer programs already prepared to 
do the scheduling. These programs vary considerably in the degree 
to which they allow consideration of job expediting costs and manpower 
availability. One class of programs, called PERT (Program Evaluation 
and Review Technique), first devised by the U.S. Navy to control 
missile design projects, go so far as to allow for probability distri- 
butions in the time estimates for individual jobs. 

The type and quality of the reports to management made available by 
critical path scheduling have caused it to be widely acclaimed as the 
best method to date for control of large, many-faceted projects. 
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4. Assembly Line Balancing: Broadly speaking, manufacturing shops can 
be divided into job shops and assembly line shops. The next topic 
indicates a use of a computer technique in job shops. Assembly line 
balancing addresses itself to the problem of work assignments along an 
assembly line so as to give personnel equal task loads and thereby 
make the overall running of the line most productive. The method also 
determines the number of workers required for a given production rate. 

Production line balancing can be of use in setting up a new line, in 
making assignment shifts for handling model changes on a set line, and 
again for shifting assignments and determining the additional workers 
needed for production rate increases and vice versa. 

The need for a computer in this technique rests on the fact that for a 
sizable assembly line the number of combinations which may arise in 
considering possible work assignments is tremendous. A definite 
sequence exists for many of the jobs to be performed. The possibilities 
of assigning these jobs, one or several to each worker, gives a large 
number of possible combinations. Added to this combinational problem 
are those jobs which can be performed any place along the line, or 
along some portion of the line's length. Finally, restrictions on such 
things as tool locations, use of the worker's left or right hand, and left 
or right side of the line must be considered. 

Mathematical techniques have been proposed and even demonstrated 
to show that an optimum set of assignments can be made. However, 
the cost of calculation for these purely mathematical approaches has 
been prohibitive. The existing successful programs for assembly 
line balancing represent a compromise between the purely analytic 
reasoning and heuristic reasoning. As such, line balancing represents 
an application of other techniques more than standing as a new technique 
by itself. 

5. Simulation: Here a mathematical and logical model of the system in 
question is created. This model is operated according to the normal 
rules of the system operation for a specified period of time. Con- 
cise data on the results is made available. This gives a picture of 
how the system will operate — quickly, and in clear, concise form. 
More important, the system may be altered or the rules changed and 
the model run again; thus, the effect of alternative decisions can be 
tested on paper. 

The term "simulation" was introduced earlier in the discussion of 
vehicle simulation. In the present context, however, the system 
simulated can be a company organization and/or its manufacturing 
tools — a much larger concept. 

Some discussion of the techniques of simulation is in order here, since, 
while the problems tackled via simulation may be extremely complex, the 
basic mathematical elements are quite simple. 

There are two basic elements of any simulation, the first being the repro- 
duction of actual (or perhaps expected) events . This can be done either 
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by using actual historical data (for example, past sales record by month) 
or by using statistical parameters which describe the desired behavior. 

To reproduce arrival of orders by size when the past 



Frequency 
Distribution Curve 




frequency distribution of order sizes is known, an accumulative distri- 
bution curve is constructed as the first step. 



100 



% of orders 
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Accumulative 
Distribution Curve 




Using random numbers (generated or taken from a table) to represent the 
"% of orders greater in size," the corresponding order size can be read 
from the curve or obtained from a table or equation representing the curve, 

This use of random numbers and the accumulative distribution curve will 
reproduce the behavior pattern from which the curve was constructed. 
Further, each order size will occur at random within the overall pattern. 
This is exactly what is desired in simulation. 

The second element in simulation of a system is time. Two distinct 
approaches have been used here: 
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1. The system is reviewed at definite intervals of time and all activities 
are updated at each review. A disadvantage of this approach is that 
for intervals of review, no activity may have occurred. However, 
many systems with natural periods can be simulated nicely using this 
procedure. Simulation of inventory and truck fleet operations fits this 
pattern. 

2. The system is reviewed only when something important happens. These 
important points in time may be recorded in a time status record 
(TSR) . Initially a distribution-of-events curve is consulted for the 
entry to the TSR. Subsequently each time an entry to the TSR is han- 
dled by updating the activity, the first step is to sample the events, 
distribution for the next TSR entry. 



Job Shop Simulation 



A job shop is distinguished from an assembly-line shop by its flexibility 
in processing articles of manufacture. The amount of time required for 
a given operation varies from article to article and the route of operations 
varies from article to article. 

Important advantages may often be realized through simulation of a job 
shop on a computer, considering such things as machines available, labor 
available, operating costs , fluctuation in orders, priority of orders, 
operating rules, etc. Records of order completions, in-process inven- 
tory, and waiting lines in front of machines can be maintained. The effect 
of proposed changes in the physical plant or operating rules may be pre- 
tested using such a simulator. 

A short description of the logic which goes into the simulation of a single 
operation job shop will illustrate the basics of this type of simulator. It 
illustrates the use of probability distributions and the time status record 
approach of handling the time element. 

The logic diagram for a single operation job shop might be as follows: 



Wait \ /Operation' 
►{ 



W 
Pieces Arrive 






The information required to simulate this system would be: 

1. The probability distribution of time between the arrival of pieces. 

2. The probability distribution of operation time per piece. 
The results desired would be: 

1. The number of pieces which can be completed per unit time. 

2. The amount of idle time experienced per unit time. 
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3. The average waiting time per piece. 

On the basis of the two probability distributions, a time status record is 
generated which guides the simulation program. The procedure for ma- 
nipulating the model is shown in Figure 15. 
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Figure 15. Job Shop Simulation 
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Two curves which might represent the required distributions are shown 
in Figure 16 along with a short, random sequence of numbers. Figure 17 
shows the time status record, and the associated recordkeeping which 
might be performed by the computer in simulating the operation. 

The conditions for beginning this simulation are: 

(a) Line 1 represents the status of the one operation shop at the beginning 
of the day. 

(b) Lines 2 and 3 are entries to the TSR made during the previous day. 

(c) Lines 4 and beyond represent entries generated during this day's 
simulation. 

(d) The simulation of this day's events begins at 9:00 a.m. (the first event) 
and can be followed in the logic chart beginning at the point marked 
START. 

This simple statistical technique for accomplishing simulation can be used 
for a wide variety of operations research problems. 
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Figure 16 
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TIME STATUS RECORD 





TIME 


EVENT 


W A F 


Zw 


Random table 
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4 
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7 
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Inventory Control 



Figure 17 



Inventory control is important to industry for sizable dollar reasons . It 
is an important area of effort computer-wise because (1) much of the 
recordkeeping required for inventory control is presently being done by 
data processing systems, and (2) many problems of inventory control can 
be tackled by those mathematical treatments which are greatly facilitated 
by use of computers. 

Inventory control, from a data processing viewpoint, is best broken into 
two distinct procedures: 

1. Forecasting the future sale of items from inventory. 

2. Use of the resulting forecast figures in a control procedure which will 
maintain that level of inventory which best meets the requirements laid 
down by management operating decisions. 

To take a hypothetical company in order to illustrate a possible study of 
these two areas, assume that the Best Parts Company has a central ware- 
house which supplies retail dealers . The situation contains the following 
factors: 

1. The inventory status must be reviewed each week. Total sales for each 
item are recorded and needed replenishments of stock are ordered at 
each review. 

2. In practice, all replenishment orders to the plant can be filled within 
a lead time of 1 1/2 weeks. 

3. The activity of different parts varies greatly. 

4. There are no largely seasonal parts. 
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5. Two years of sales records are available. This record gives the 
total sales of each part for each one-week period. 

Best Parts Company management believes that it can improve the control 
of the inventory in one or more of three ways: (1) reducing inventory, 
(2) reducing out-of-stock conditions, and (3) giving simpler processing 
of control information. 

A study of the control procedure is initiated with the understanding that 
the computer facility, with whatever new techniques available, may be 
used. The study takes the following form: 

Forecasting 

Four approaches to arriving at the forecast sales to retailers for each 
part are considered: 

1. The use of linear regression in finding a prediction relationship of the 
form 

Forecast Sale = a + a^x^ + a£X2 + a n x n 

where the x's may be factors felt to correlate with or influence the 
sales, such as gross national product, number of people in $5,000 - 
$10,000 income bracket, etc. This approach requires extensive 
analytic effort, much diverse data, and large amounts of computing 
effort. 

2. Classical time series analysis gives a forecast based only on analysis 
of previous sales,data. The linear trend line of the form 

Forecast Sale = a + bT 

where T is time in sales periods, requires only past sales data, but 
considerable computing for fitting of a and b. 

3. A common practice in industry is to use a moving average for sales 
forecasting. The formula for a four-week moving average is 

F ., = S + S , + S „+S n 
r+1 r r-1 r-2 r-3 



where F , is the forecast sales for week r+1 and the S's are the 
actual sales for the previous four weeks. This method is simple, 
data-wise and computation- wise . A major disadvantage is that it does 
not generate statistics which are found important in the control pro- 
cedure — in particular, the statistic called the standard deviation of 
actual sales from the forecast. Approaches 1 and 2 may generate 
this statistic. 

4. A weighting technique credited to Robert K. Brown and called exponen- 
tial smoothing has gained attention and application in the past few years . 
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Forecasting is based, as in the moving average, on past sales data 
and requires even less computing. A major advantage is that a sim- 
ply computed estimate of the required standard deviation is available 
and corrected to the most recent sales review figures. 

The exponential smoothing procedure is chosen for study on the basis of 
limited available data and limited computing time . The forecast equations 
for exponential smoothing may be given as 

Pj = P^ +*T (S i _ 1 - Pi*_i) (Demand) 

T t = cC (Pj - Pi*_!) + (1 - oC ) T._ 1 (Trend) 

Pj* = Pj + ® sZ ^ (ly (Forecast) 

D i = * < S i-l " P i*-1> + d - °< ) D i-i (Deviation) 
where 

S- -, = sales figure for period just ending. 



Then 



and 



with 



P^ = demand expected for next period 

Tj = trend measure of whether last several periods have seen a 
rise or fall 



P|* = trend adjusted, final forecast, for next period sales 



1. 25 Dj = estimate of standard deviation of actual sales from fore- 
cast. 

The constant o< represents the weighting choice . Without the trend ad- 
justment, the demand (Pi) computed will closely approximate the moving 
average based on N periods if oC is chosen as 

2 

With the use of the trend adjustment, experience has shown that a choice 
of cK = . 1 is most satisfactory. 

Control Procedure 

The key to inventory control is in the way one parameter is computed. 
This parameter is the requirements figure. The requirements figure 
must be arrived at in such a way that it relects (1) the knowledge of actual 
sales, so that stock is available to meet the demand, and (2) management's 
ideas of what is adequate by computing the proper stock level, within a 
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►J 



Forecast 



Lead 



defined approximation, so that out-of-stock conditions are kept below a 
set level. 

This implies that management must set the level for out-of-stock condi- 
tions — and this is true. The familiar quality decisions such as "we must 
have few out-of-stock conditions" or "we must reduce inventory" must be 
replaced with decisions such as "an average of 5% out-of-stock is to be 
our goal." 

Literature is available on methods of arriving at a proper out-of-stock 
percentage figure which tends to minimize inventory operating costs for 
particular inventory situations. 

In this example, the management percentage level figure is assumed 
available . 

One approach to defining the necessary requirements equation is to con- 
struct an inventory model similar to that shown below. This model assumes 
the following "ideal" conditions: 

1. It represents a possible inventory level as a function of time for one 
item. 

2. The sales are linear — the same sales figure for each period. 

3. The sales in each period are exactly the amount forecast for that 
period. 

4. The lead time is 1 1/2 review periods (lead time being the time be- 
tween generation of an order and actual receipt of replenishment 
stock) . 

The graph begins at a point in time where the actual inventory is 
exactly at the requirements level. 



Requirements Level 




Time 
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The upper curve represents the inventory actually in the warehouse 
at any given time plus the inventory on order. 

The graph shows that for a requirements figure computed strictly in terms 
of the forecast and the lead time, zero stock conditions are to be expected. 
Since the Best Parts Company knows that it is not operating under "ideal" 
conditions (sales fluctuate from period to period), a safety stock is needed. 
The final requirements equation then is: 

R i = p i* + Lead * p i* + Safety 
The Safety Stock Calculation 

Reference must be made here to the discussion of probabilities in the 
chapter on empirical relationships. If a frequency distribution of errors 
in the forecast over a period of time is plotted, it might look like Figure 
18 (assuming a normal distribution) . 

The accumulative area under the curve, as we move from left to right 
along the error axis , represents the probability of the error in forecast 
being equal to or less than the corresponding value of the error. 




-2-10 l 

Error 
(Forecast — Actual Sale) 

Figure 18 

If the shaded area in Figure 18 represents 5% of the total area, then" only 
5% of the time is the error in forecast greater than four units. 

Therefore, if four units are always added to the inventory stock require- 
ments , the out-of-stock conditions should occur only for errors larger 
than four units — or only 5% of the time . 

For computational purposes, the same proper error value for a stated 
probability may be obtained by using the transformation equation: 



t = 



_ x-m 



where for this use 
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t = value corresponding to required probability level and is taken 
from the tabulated values available for the normal distribution 

x = desired error value 

cr = standard deviation of error 

m = mean of the errors (in this case m = o) 

Then 

x = cr t 

and since 1. 25 Dj_-|_ is an estimate of cr, the safety stock required may 
be computed as 

x = 1.25 Dj_ *t 

The value of t corresponding to a probability level allowing 5% out-of-stocks 
is 1.65. The final requirements equation to be used by Best Parts Com- 
pany is 

Rj = P.* + Lead • Pj* + 1.65 • 1.25 • D i 

Note that the statistic D, like all the other quantities here, is computed 
for each item in the inventory. The requirements figure for an item level 
is thus dependent on that item's own past sales behavior. Items with steady 
sales will have small safety stocks. Items with widely fluctuating sales will 
have large safety stocks. This may often make possible the seeming para- 
dox of reducing overall inventory and at the same time reducing out-of-stock 
situations . 

This completes the two basic selections needed in establishing an inven- 
tory control procedure — a forecasting technique and a requirements com- 
puting technique. All other quantities relating to inventory control can be 
computed each period in terms of the forecast and the requirements. For 
example, the replenishments equation for the inventory might be given as: 

ORDER = REQUIREMENTS - ON HAND - ON ORDER 

This procedure requires extensive testing before it may be successfully 
initiated. One way of accomplishing this testing is to simulate on a com- 
puter what would have happened if the procedure had been used over some 
time interval in the past. A simulation of inventory control has two basic 
elements: actual events and time. Actual sales data (quantity of sales 
of each item or each period in the time interval) is all that is needed to 
satisfy the reproduction of actual events elements . The normal use of a 
review period in inventory control makes it a natural for the use of re- 
view at intervals in handling the time element. 

A logic flow chart for a computer program to simulate the inventory model 
for a single inventory item is shown in Figure 19. 



145 





KEY 


Rq 


= Receipts for Period i 


ORi 


= On Order Total 




Beginning Period i 


OHj 


= On Hand Beginning 




Period i 


N 


= Number of Periods 




to Simulate 


L- 


= Lead Time in 




Periods 



Set Starting Values for 

p i*-l» Si.!, Tj.!, D.^ 
OR. ., OH. . 

RCi_i, RC., Rc i+L-1 



COMPUTE AND 
RECORD 
OH^OHj.!- 
Sj-i + RCj.! 




COMPUTE AND 
RECORD 

Pi = P i *.l+«<(Sn-P i *. 1 ) 
Ti =^<(P i -P i *_ 1 )+(l-K)T i _ 1 

p i* = p 1 + (i 5r ,T i 

Di^Si-P^jHlwODi.! 
Ri = Pj* + L(Pj*) + 2.45Dj 




COMPUTE AND 
RECORD 

R C i+ L = Ri-°Hi + 
ORi-1 



COMPUTE AND 
RECORD 
OR. = OR. 




RECORD 
OUT OF STOCK 



Set to Next 

Period 
i = i + 1 



Set 
RC i + L = 



END 



Figure 19 
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An important result of such a simulation is the data needed to prepare 
charts similar to that shown in Figure 20 . 



800 _ 



600 _ 



Number 

of 
Pieces 



400 



200 - 



X Period Sales 

On Hand 

A Requirements = On Hand + On Order 



,A. 



&—--£.. 



.^. 



,A— — fr- 



-•&—-. 



*' 



'^r 



„— A 



'-*' 



'-*' 




Periods Ending 
Figure 20 



A graphic display of this kind, showing what will happen if a procedure is 
adopted, represents an excellent pretesting of the practicality of the pro- 
cedure . 



Linear Programming 



Mathematically, linear programming is finding a solution of a set of 
equations similar to: 

a ll x l + a 12 x 2 + a 13 x 3 + a 14 x 4 + a 15 x 5 = b l 
a 21 x l + a 22 x 2 + a 23 x 3 + a 24 x 4 + a 25 x 5 = b 2 

a 31 X l + a 32 X 2 + a 33 X 3 + a 34 X 4 + a 35 X 5 = b 3 

in which we have more variables than equations, and therefore an infinite 
number of solutions. 

The particular solution chosen is the one which maximizes or minimizes 
an additional equation (generally a cost equation) which might look like: 

cost = c.x-, + c 2 Xo + CgXo + C A X 4 + C 5 X 5 

Any system (either physical, such as a petroleum refinery, or operational, 
such as work assignments, or a combination of both) which may be de- 
scribed in terms of linear relationships may be studied by linear program- 
ming. 

Many relationships which are not linear may be approximated by combi- 
nations of linear segments as illustrated graphically in the accompanying 
figure . 
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f(x) 




Much of the present-day power of linear programming is due to its having 
been taken out of the hands of the mathematicians and reduced to a fairly 
common engineering procedure. Rather than working with sets of equa- 
tions, the engineer works with an array or matrix which is the format re- 
quired for entering data to the computer. In place of mathematical opera- 
tions in creating this matrix, the engineer employs straightforward rules. 



.Variables 
Restrictions ^v_ 


Straight 

Run 

Gasoline 


Alkylate 


Base 
Stock 


Butane 


Re form ate 
One 


Reform ate 
Two 


Crude 
One 


Crude 
Two 


Lead 
One 


Lead 
Two 


Lead 
Three 


Excess 


c 
S 

•3 

a" 

OS 


Profit Coeff. 


S.O 


5.0 


5.0 


5.0 


4.80 


4.60 


-2.80 


-3.10 


-.092 


-.092 


-.092 


-10.0 




1 Octane 


.3 


-.7 


-.1 


-.09 


-.1 


-.3 






-.1 


-.3 


-.7 







2 Boiling Point 


.8 


.3 




.248 


-.7 


-.7 

















3 Vapor Pressure 


1.0 


-.5 


-1.0 


2.6 


-2.0 


-2.0 

















4 Quantity 


1.0 


1.0 


1.0 


.1 


1.0 


1.0 












-1.0 


1.0 


S Lead 1 


-l.S 


-1.5 


-1.5 


-.15 


-1.5 


-1.5 






1.0 











6 Lead 2 


-1.0 


-1.0 


-1.0 


-.1 


-1.0 


-1.0 








1.0 









7 Lead 3 


-.S 


-.5 


-.5 


-.05 


-.5 


-.5 










1.0 







8 Reformer Charge 










1.11 


1.25 


-.15 


-.2 










.5 


9 Reformer Capacity 










1.11 


1.25 


.05 












.6 


10 Butane 








.1 


-.03 


-.04 




-.01 










.05 


11 Straight Run 
Gasoline 


1.0 












-.1 


-.15 










.15 



Figure 21. 



A look at the array in Figure 21, which represents a complex problem to 
be optimized, will illustrate the use of such rules. 

This array represents the statement of the problem to be optimized. The 
objective will be to determine a recipe for blending available petroleum 
components to obtain a motor fuel meeting several specifications — and 
return maximum profit to the refinery. 

Each horizontal row in the array represents a restriction (equation) — and 
each vertical column a variable. The actual numbers in the array are 
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coefficients which would appear in the corresponding set of linear equa- 
tions which are required to describe the problem mathematically. The 
variables are the components to be blended. The restrictions may be of 
two types: 

1. Quantity restrictions — for example, where the amount of motor fuel 
to be produced is restricted. 

2. Qualitative restrictions — for example, where the Reid vapor pressure 
of the resulting fuel is restricted. 

Line 4 in the body of the array is read: 1.0 times amount of straight run 
gasoline to be used in the blend, plus 1.0 times amount of alkylate to be 

used in the blend, plus , minus 1.0 times excess motor fuel, 

must be equal to the 1. barrels of motor fuel which is required. The 
coefficients can be entered to the array without ever putting down the 
equation stated above for this quantitative restriction. 

Row 1 represents the representation of the octane number restriction. 
The refinery sets a minimum level for the resultant octane number. In 
this case, let the minimum specification be S and let the actual octane 

ratings of the variables X^ to be blended be 1^ for i = 1, 2 n 

where n is the number of variables . 

The assumption is that the resultant blend octane rating will be a weighted 
average by volume of the individual specification or: 

n 
S + P = Z Ii Xi 
i=l 



n 




n 


z 


(S - Ii) Xi + P 


i 


i=l 




i=l 



n 

Z Xi 
i=l 

where the P allows for the weighted average to be greater than the mini- 
mum specifications. A suitable form of the above equation is: 



Xi = 



n 
The term P Z Xj is called the slack. 
i=l 

Note that the coefficients for each of the variables is to be S-Ii or the 
minimum specification minus the individual specification. 

This characteristic allows two rules for automatically entering coefficients 
into an array to handle most quality restrictions: 

1. The magnitude of the coefficient is to be the absolute value of the 

difference between the required specification and the individual speci- 
fications . 
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2. The sign of the coefficient is chosen positive if the individual specifi- 
cation exceeds the required specification, negative if the individual 
specification is below the required specification. 

Linear programming codes exist for most of the digital computers in 
commercial use. Once the array has been determined to adequately de- 
scribe the system under study, the data represented by the array is pre- 
sented to a computer along with the linear program . 

The output to be expected from the optimization computation is: 

1. Variable activities. This is a list of those variables chosen to satisfy 
the restrictions . With each variable is given the amount of the variable 
to be used in the resultant product. 

Many linear programs also indicate the range, if any, through which 
this amount may be varied, leaving the solution at a maximum or 
minimum as the case may be. 

2. Shadow prices. For those variables not selected as part of the solu- 
tion, an indication of the reason is given by the shadow prices. The 
shadow price for an unselected variable represents the cost of forcing 
one unit of the variable into the solution, considering the adjustments 
that would be necessary on other variables to again satisfy the stated 
restrictions. Again there is often indicated a range through which 
the cost or profit coefficient for the unselected variable may be moved 
without removing the optimization. 

Results of linear programming computation are used as a recipe for 
day-to-day selection of alternatives in such diverse areas as gasoline 
blending of available stocks, mixture of grain for feed, blending of 
meats for sausage, and slitting of paper rolls and fabric stock. 

The technique is also used on a study basis to aid long-range planning 
in the selection of refinery crude stocks, the selection of machine 
components in process industries, the scheduling of machine loads, 
and other areas . 

SUMMARY ON OPERATIONS RESEARCH 

As a brief summary of Operations Research, there are several techniques 
which have gained considerable acceptance and use in the study of man/ 
machine systems: linear programming, forecasting, and simulation. 

Linear programming computes the optimum choice of decision alterna- 
tives for a linear system within the available definition of that system. 

Forecasting by statistical means has important uses in long-range plan- 
ning and in day-to-day control of such things as inventory. 
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Exercise 



Simulation may be used to study systems thus far undescribed mathemat- 
ically, by reproducing the logical operation of the system in testing alter- 
nate decisions. 

Critical path technique for scheduling is used to control and improve the 
planning, scheduling, and expediting of projects . 



1 . A company would like to simulate the operation of a fleet of trucks 
which deliver from one factory to many parts of the country. The 
delivery is to be made for one product. The company has the follow- 
ing data which is important to the simulation: 



Day 



Number of Orders 
(Units of 1 Truckload) 



100 
97 

115 
40 
80 



150 



102 



Number of Days 
Required to Deliver 
and Return 



Number of Times 

This Period Was 

Required 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 



200 

500 

700 

2000 

4000 

5500 

3500 

3000 

500 

300 

40 





a. Sketch the distribution curve which might be used in simulation of 
the receipt of orders. 

b. Sketch the distribution curve which might be used to determine the 
length of time required to deliver each order (turnaround time) 
during the simulation. 
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The company would like to determine the best number of trucks to use 
to satisfy the delivery requirements while keeping the truck fleet cost 
as low as possible. 

The company estimates that it costs D dollars to maintain and operate 
one truck each year, and that the cost of back orders in terms of future 
sales lost is: 

C = AN X + BN 2 

where A and B are constants, N^ is the number of back orders not 
more than seven days old, and N 2 is the number of back orders over 
seven days . 

Deliveries may be made every day of the year . 

c. Draw a logic flow chart of a simulation program to operate for one 
year, which the company might use to determine the optimum num- 
ber of trucks to maintain in the truck fleet. 
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EXERCISE 6 



Generalized flow chart to compute F(x) = a^ x n + a2X n_1 +a. 



n+1 





READ 










N, A. 






1 = 1,2- 


— ,N+1 






V V 






READ 






X 






V 






Set 






1= 1 






F(x) 


= 0. 
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r i 


r 










Compute 




Step 




F(x) = F(x)«X + A. 




1 = 1+ 1 




i 


i 




1 


' 






^^Test ^s. /^~X 




S^I = N+ 1 ?^^ iiiu >■ 




r 


1 


r 




PRINT 






X, F (x) 
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Section IV 



1 C FOR COMMENT 




ISUUMJNT * 
* NUMttR ,_ 

1 3 t 


FORTRAN STATEMENT 

7 10 15 30 23 30 35 *0 «3 30 33 «0 63 70 72 


c. . . . 


| B i-Efi-xC 1 S£ ,// , , , , , 






I'.l , I I I ................. 






x = . o,r 






SUMly- O . 0'--\ I I I 1 I I . 






SUM<A*- O .0 , , . , , , , , , 




. . . .1 


Y*l.. i O/.<l.. t D+X*)t i ) ......... . .... 






6.0. t,4 .(«„2.). t .f. , . . , . . , , ............. 




... .4 


suM*J s = suM4-ty. . . , . , , , , , , ........ 






1.--.2 ............................. .......... 






6.0. .T,0 1 , 




. . . 1 


SUMZ s -SUMZ^-f-y, . . , , , i , , 






l.\l. .......................... . . . . . 




3 


X=.X>,.<?.r. ...... ........ ....... 






lF.(X-l.0\l,Si.^ i ............... . ...... 




. S 


sUM-,0. OS/,3. 0-*.(,l.S*4. 0*SU,M^+2.,0*SUM i 2) 






PRI^TZO^SJJM I,,,,, 






£X^ ..................... . . 
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Section V 



EXERCISE 3 

The initial conditions may be used to complete the solution at t = . 



<*0 = 



dt 



B --v(£) 



(S) o - X o^o(£) o 



Succeeding values for the variables and their derivatives may be obtained 
by two-level iteration. For each value of t, iteration is used to find a 
solution to o< = t + sin a(. t. The following logic diagram makes this clear: 



Complete 
initial solution 



Step 
'n+^'n+At 



i__£ 



= 0, 1, 2, 3- 





F(K) =«| - t - linKt 
F' (*j) = 1 - eta «, t n 



H = «i+l 

1 



= 0, 1, 2, 3, 4 





Set 






"» S «M 






1 




Compute 




©.-<.?l-<§L" 


(Sl-»-„-(£1 


(§) n --vn-„(^) n 
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Section VI 






EXERCISE 2 




Inverse is 




-4 4 -1 




4-5 2 




-1 2 -1 




EXERCISE 3 




x 1 = 2 




x 2 = 5 




x 3 = 9 




EXERCISE 4 



Solution is 
x, = 3 



x 2 =-3 



/28 


/28 


/28 


/31 


/37 


50 

41 


33 


34 


29 


30 


29 


28 


28 


28 


26 


28 


28 


28 


27 


25 
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Section VII 

EXERCISE 3 



t = x - m = . 0027 - .0025 = 1>0 
°" .0002 



for t = 1.0, from table of probabilities 

00 

2 f -t 2 /2 
' ^ e dt = . 1587 



27T 



2.5 
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APPENDIX II : BIBLIOGRAPHY 



The Data Processing Bibliography (form J20-8014) provides abstracts of 
several excellent texts on numerical techniques, scientific applications 
and related subjects in the field of data processing. IBM publications of 
interest to the engineer include: 

Programmer's Primer for FORTRAN Automatic Coding System for the 
IBM 704 Data Processing System (F28-6019) 

Random Number Generation and Testing (C20-8011) 

Flow Charting and Block Diagramming Techniques (C20-8008) 

Introduction to IBM Data Processing Systems (F22-6517) 

Inventory Management Simulation (E20-8063) 

AUTOPROMT — Automatic Programming of Machine Tools (E 20-6 104) 

Production Line Balancing at Westinghouse (E20-4037) 

PERT — A Dynamic Project Planning and Control Method (E20-8067) 

Optimum Shop Scheduling and Machine Loading (E 20-4044) 

Quality Control (320-0955) 

Improved Job Shop Management through Data Processing (E20-6071) 
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